Typical and rare fluctuations in nonlinear driven diffusive systems with dissipation 
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We consider fluctuations of the dissipated energy in nonlinear driven diffusive systems subject to 
bulk dissipation and boundary driving. With this aim, we extend the recently-introduced macro- 
scopic fluctuation theory to nonlinear driven dissipative media, starting from the fluctuating hydro- 
dynamic equations describing the system mesoscopic evolution. Interestingly, the action associated 
to a path in mesoscopic phase-space, from which large-deviation functions for macroscopic observ- 
ables can be derived, has the same simple form as in non-dissipative systems. This is a consequence 
of the quasi-elasticity of microscopic dynamics, required in order to have a nontrivial competition 
between diffusion and dissipation at the mesoscale. Euler-Lagrange equations for the optimal density 
and current fields that sustain an arbitrary dissipation fluctuation are also derived. A perturbative 
solution thereof shows that the probability distribution of small fluctuations is always gaussian, as 
expected from the central limit theorem. On the other hand, strong separation from the gaussian 
behavior is observed for large fluctuations, with a distribution which shows no negative branch, 
thus violating the Gallavotti-Cohen fluctuation theorem as expected from the irreversibility of the 
dynamics. The dissipation large-deviation function exhibits simple and general scaling forms for 
weakly and strongly dissipative systems, with large fluctuations favored in the former case but 
heavily supressed in the latter. We apply our results to a general class of diffusive lattice models 
for which dissipation, nonlinear diffusion and driving are the key ingredients. The theoretical pre- 
dictions are compared to extensive numerical simulations of the microscopic models, and excellent 
agreement is found. Interestingly, the large-deviation function is in some cases non-convex beyond 
some dissipation. These results show that a suitable generalization of macroscopic fluctuation theory 
is capable of describing in detail the fluctuating behavior of nonlinear driven dissipative media. 



I. INTRODUCTION 

Fluctuations are inherent to many physical phenom- 
ena, reflecting the hectic microscopic dynamics at macro- 
scopic scales. In spite of their apparent random ori- 
gin, essential physical information is encoded therein pQ. 
A classical example is the fluctuation-dissipation theo- 
rem, which relates the linear response of a system to 
an external perturbation to the fluctuation properties 
of the system in thermal equilibrium [21 [3]- More re- 
cently, the investigation of general properties of fluctu- 
ations in nonequilibrium steady states is opening new 
paths for understanding physics far from equilibrium j3J- 
H4] , The study of fluctuation statistics of macroscopic 
observables provides an alternative path to obtain ther- 
modynamic potentials, a complementary approach to the 
usual ensemble description. This observation, valid both 
in equilibrium p] and nonequilibrium 19] , is most rele- 
vant in the latter case because no general bottom-up ap- 
proach, connecting microscopic dynamics to macroscopic 
nonequilibrium properties, has been found yet. In this 
way, the large deviation function (LDF) controlling the 
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statistics of these fluctuations may play in nonequilib- 
rium statistical mechanics a role similar to the equilib- 
rium free energy [H [TU] . A central point for this emerg- 
ing paradigm is the identification of the relevant macro- 
scopic observables characterizing the out of equilibrium 
behavior of the system at hand. The system dynamics 
often conserves locally some magnitude (a density of par- 
ticles, energy, momentum, charge, etc.), and the essential 
nonequilibrium observable is thus the current or flux sus- 
tained by the system when subject to boundary-induced 
gradients or external fields. Therefore, the understanding 
of current statistics in terms of the microscopic dynam- 
ics represents one of the main problems of nonequilib- 
rium statistical mechanics, triggering an intense research 
effort in recent years. In this context, key results are 
the Gallavotti-Cohen fluctuation theorem [5], which re- 
lates the probability of observing a given current fluctu- 
ation J with the probability of the reversed event — J, or 
the recently-introduced Isometric Fluctuation Relation 
[14] , which relates the probability of any pair of isomet- 
ric current fluctuations (J, J'), with |J| = \J'\. These 
intriguing "symmetries" appear as a consequence of the 
invariance under time reversal of the underlying micro- 
scopic dynamics. These and other recent results [5rTB] 
are however restricted to nonequilibrium "conservative" 
systems characterized by currents. 

On the other hand, many nonequilibrium systems are 
inherently dissipative, that is, they need a continuous in- 
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put of energy in order to reach a steady state. In this class 
of systems, the relevant macroscopic observable is not 
only the current: the dissipated energy is also expected to 
play a main role. These systems include granular media 
[H] , dissipative biophysical systems [TH] , turbulent fluids 
[19], active matter [20], chemical reactions [21], popula- 
tion dynamics |22j . etc. In general, this class comprises 
all sort of reaction-diffusion systems where dissipation, 
diffusion and driving are the main physical mechanisms. 
Fluctuations in dissipative media have been much less in- 
vestigated, most probably because their physics is more 
complicated as a result of the irreversibility of their mi- 
croscopic dynamics. In principle, most of the results for 
nonequilibrium steady states we have referred to in the 
previous paragraph are not applicable to dissipative sys- 
tems, since they stem from the reversibility of the un- 
derlying microscopic dynamics in the conservative case. 
Therefore, a question naturally arises as to whether it is 
possible to extend some of these ideas to dissipative me- 
dia. One of the main goals of the present paper is to give 
a (partial) answer to this question. 

In this work, we analyze both typical and rare fluc- 
tuations in nonlinear driven diffusive systems with dis- 
sipation. This is done by combining a suitable general- 
ization of macroscopic fluctuation theory (MFT) [7] to 
the realm of dissipative media, and extensive numerical 
simulations of a particular albeit broad class of micro- 
scopic models. Our starting point is a general fluctu- 
ating balance equation for the (energy) density, with a 
drift term proportional to the spatial derivative of the 
current and a sink term. This mesoscopic description is 
expected to be valid for many driven dissipative media 
over a certain "hydrodynamic" time scale, much larger 
than the one characteristic of the microscopic dynam- 
ics. Over the fast (microscopic) time scale, the system 
forgets the initial conditions and relaxes to a local equi- 
librium state in which all the properties of the system 
become functionals of a few "hydrodynamic" fields, here 
the density, the current and the dissipation. Afterwards, 
over the much slower hydrodynamic time scale, the sys- 
tem eventually approaches the steady state following the 
mesoscopic balance equation. We focus on the fluctua- 
tions of the system in this nonequilibrium steady state, 
in which the dissipation and the injection of energy bal- 
ance each other. By using this fluctuating hydrodynamic 
picture together with a path integral formulation, we de- 
rive a general form for the action associated to a history 
of the density, current and dissipation fields (that is, a 
path in mesoscopic phase space) . Remarkably, this action 
takes the same form as in conservative nonequilibrium 
systems [THE] , simplifying the analysis in the dissipative 
case. This is both an important and a surprising result, 
which stems from the quasi-elastic character of the un- 
derlying microscopic dynamics in the large system size 
limit. This quasi-elasticity is necessary in order to have 
a balanced competition between diffusion and dissipation 
at the mesoscopic level. From the derived action func- 
tional, and using the recently-introduced additivity con- 



jecture [51 ITTI - ITS] . a general form for the LDF of the dis- 
sipated energy is derived, with a "Lagrangian" including 
second order derivatives. Therefrom, we derive the Eulcr- 
Lagrange equation (a fourth-order differential equation) 
for the optimal fields responsible of an arbitrary fluc- 
tuation. This Lagrangian variational problem can be 
mapped onto an equivalent Hamiltonian problem (four 
coupled first-order differential equations) which turns out 
to simplify the analysis. We use this Hamiltonian pic- 
ture to analyze in detail three different limits, namely 
small fluctuations around the average for arbitrary dissi- 
pation coefficient, and the whole spectrum of fluctuations 
(typical and rare) for weakly- and strongly-dissipative 
systems. The statistics of typical (that is, small) fluc- 
tuations is gaussian as expected from the central limit 
theorem. However, strong separation from gaussian be- 
havior is observed for rare fluctuations, with a distri- 
bution which shows no negative branch, thus violating 
the Gallavotti-Cohen fluctuation theorem as otherwise 
expected from the irreversible character of microscopic 
dynamics. We study in general the weakly-dissipative 
system limit using a singular perturbation expansion. 
This yields a simple scaling form for the dissipation LDF, 
showing that large dissipation fluctuations are favored in 
this weakly-dissipative limit, with a LDF which extends 
over a broad regime and decays slowly in the far positive 
tail. On the other hand, a different perturbative analy- 
sis in the strongly-dissipative system limit can be carried 
based on the formation of boundary energy layers in this 
limit for all fluctuations, which effectively decouples the 
system in two almost-independent parts. This analysis 
shows that large dissipation fluctuations are heavily sup- 
pressed in this limit, as opposed to the weakly-dissipative 
regime result. 

We apply this theoretical scheme to a general class 
of d-dimensional dissipative lattice models with stochas- 
tic microscopic dynamics, for which the hydrodynamic 
fluctuating picture used above as starting point can be 
demonstrated in the large system size limit |39j (we will 
focus here in one dimension for simplicity). In these mod- 
els there is one particle at each lattice site, characterized 
by its energy. Dynamics is stochastic and proceeds via 
collisions between nearest neighbors, at a rate which de- 
pends on the energy of the colliding pair. In a collision, 
a certain fraction of the pair energy is dissipated, and 
the remaining energy is randomly distributed within the 
pair. This mechanism gives rise to a nonlinear competi- 
tion between diffusion and dissipation in the macroscopic 
limit, provided that the microscopic dissipation coeffi- 
cient scales adequately with the system size. This class of 
models represents at a coarse-grained level the physics of 
many reaction-diffusion systems of technological as well 
as theoretical interest. In particular, when the collid- 
ing pair is chosen completely at random, independently 
of the value of its energy, the Kipnis-Marchioro-Presutti 
(KMP) model [23] for heat conduction is recovered in the 
conservative case. The KMP model plays a main role in 
nonequilibrium statistical physics as a touchstone to test 
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theoretical advances [5J [51 [HI ITTHT51 123"] . Our general 
class of models contains the essential ingredients char- 
acterizing most dissipative media, namely: (i) diffusive 
dynamics, (ii) bulk dissipation, and (hi) boundary injec- 
tion. The chances are that our results remain valid for 
more complex dissipative media described at the meso- 
scopic level by a similar evolution equation. Here we 
report analytical and simulation results for the statistics 
of the dissipated energy in this general class of models 
using both standard simulations and an advanced Monte 
Carlo method [27J . The latter allows the sampling of the 
tails of the distribution, and implies simulating a large 
number of clones of the system. 

The plan of the paper is as follows. Section[IT]describes 
a suitable generalization of macroscopic fluctuation the- 
ory to nonlinear driven dissipative media. The large- 
deviation statistics of the dissipated energy, a central 
observable in this type of systems, is investigated here 
within both the Lagrangian and Hamiltonian equivalent 



frameworks. Section III is devoted to the detailed study 
of different asymptotic behaviors within the Hamiltonian 
formulation, which turns out to simplify the analysis. In 
section [TV] we define a general class of microscopic lattice 
models whose stochastic dynamics is dissipative [351 EU • 
The theoretical framework developed in the previous sec- 
tions is applied to this family of models in section [Vj and 
the LDF for the dissipated energy is explicitly worked 
out. The analytical predictions are compared to exten- 
sive numerical simulations of the microscopic models, and 
a very good agreement is found. A summary of the main 
results of the paper, together with a physical discussion 
thereof, is given in sec. |VI| Finally, the appendix deals 
with some technical details that, for the sake of clarity, 
we have preferred to omit in the main text. 



II. MACROSCOPIC FLUCTUATION THEORY 
FOR DRIVEN DISSIPATIVE SYSTEMS 

In this work, we will analyze a general class of sys- 
tems whose dynamics at the mesoscale is described by 
the following fluctuating evolution equation 



d t p(x,t) = -d x j{x,t) + d(x,t) 



(2.1) 



We focus here in one dimension for simplicity, but our 
analysis can be carried out in an equivalent manner in d- 
dimensions. Ineq. (2.1 1, p(x, t), j(x t t) and d(x, t) are the 



density, current and dissipation fields, respectively, and t 
and x e [—1/2, 1/2] are the macroscopic time and space 
variables, obtained after a diffusive scaling limit such that 
x = x/L and t = t/L 2 , with x and t the microscopic 
space and time variables and L the system length. These 
coarse-grained spatial and temporal scales emerge from 
a suitable continuum limit of the underlying microscopic 
dynamics |39j . The current field is a fluctuating quantity, 
and can be written as 



The first term is Fourier's law, where D(p) is the diffu- 
sivity (which might be a nonlinear function of the local 
density), and t) is the current noise that is gaussian 
and white, 

(£(*, *)) = 0, t)ttf, 0) = ^6(x-x')6(t-l/), 

(2.3) 

with a(p) being the so-called mobility. This gaussian 
fluctuating field is expected to emerge for most situa- 
tions in the appropriate mesoscopic limit as a result of a 
central limit theorem: although microscopic interactions 
for a given model can be highly complicated, the ensuing 
fluctuations of the slow hydrodynamic fields result from 
the sum of an enormous amount of random events at 
the microscale which give rise to gaussian statistics, with 
an amplitude of the order of L" 1 ^ 2 , in the mesoscopic 
regime in which eq. (2.1) emerges. On the other hand, 



the dissipation field d(x, t) is 



d(x,t) = —vR(p(x,t)), 



(2.4) 



where v is the macroscopic dissipation coefficient, and 
R(p) is a certain function of the density p. For the calcu- 
lations which follow throughout this section, it is useful 
to introduce a new variable y, such that 



y = R(j>), 



d(x, t) = —vy(x, t). 



(2.5a) 



(2.5b) 



The dissipation field is present at the mesoscopic level 
because the microscopic stochastic dynamics of the mod- 
els of interest dissipates some energy, that is, we have 
the equivalent of a microscopic restitution coefficient a, 
so that the amount of dissipated energy is proportional 
to 1 — a. The macroscopic dissipation coefficient v is 
thus proportional to 1 — a. Note however that there is 
no noise term in eq. (2.4 1, so the local fluctuations of 



the dissipation field are enslaved to those of the density 
p(x, t). The physical reason for this behavior is that the 
microscopic dynamics must be quasi-elastic in order to 
ensure that dissipation and diffusion take place over the 
same time scale in the thermodynamic limit. Typically, 
1 — a must scale as L~ 2 that is the order of magnitude 
of the diffusive term in a system of length L [251 132] ■ 

The boundary conditions for eq (2.1 1 depend on the 
physical situation of interest. For instance, we may con- 
sider that the system is kept in contact with two ther- 
mal reservoirs at x = ±1/2, at the same temperature T, 
so /o(±l/2, t) = T. In that case, the system eventually 
reaches a steady state in the long time limit, for which 
the injection of energy through the boundaries and the 
dissipation balance each other. The stationary average 
(macroscopic) solution of (2.1) verifies 



j(x,t) = -D{p)d x p(x,t) + £(x,t). 



f av (x) + uR{p m {x)) = 0, j av {x) = -D(p av (x))p' av (x), 
where the prime indicates spatial derivative. The first 



(2.2) equation in (2.6) follows from (2.1), and the second one 
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is Fourier's law for the averages. Equivalently, a closed 
second-order equation for p may be written, 



d_ 

dx 



[D(p av )p' av } = vR{p a 



(2.7) 



with the boundary conditions p av (±l/2) = T . Equations 
(2.6) and (2.7) can be also written for the variable y 
introduced in cq. ( |2.5a l, 



j' av (x) + yy av {x) 
with 



0, Ja 



-D(y av (x))y' av (x), 
(2. 



D(y) 



since 



j av (x,t) = -D(y av )d x y(x,t). 



(2.9) 



(2.10) 



Thus D is an "effective" diffusivity: it is the factor multi- 
plying the spatial gradient when writing Fourier's equa- 
tion in terms of the new variable y. Equation (2.8) can 



also be summarized in a second order differential equa- 
tion for y av , 



D{y av )y' a , 



vy a 



y av (±l/2) = R(T). (2.11) 



Interestingly, it can be shown (see below) that D is con- 
stant, independent of y, whenever y = R(p) depends al- 
gebraically in p, a case we will study in detail in section|V] 
This observation considerably simplifies the subsequent 
analysis. 

The probability of observing a history {p(x, t),j{x, t)}^ 
of duration r for the density and current fields, starting 
from a given initial state, can be written now as a path 
integral over all the possible realizations of the current 
noise {£(x,£)}q, weighted by its gaussian measure, and 
restricted to those realizations compatible with eq. ( |2.1[ ) 
at every point of space and time |12j . This probability 
hence obeys a large deviation principle of the form [U [7J- 



P({p,jVo)~exp(+Ll T {p,3}). 



(2.12) 



with a rate functional [71 [5] 
Zt[p,3\ = - I dt 



-1/2 2cr(p) 



with p(x,t) and j(x,t) coupled via the balance equa- 
tion ( |2.1[ ), and the dissipation d(x,t) given in terms of 
p(x,t) by (|2.4[). Equation (2.13) expresses the gaussian 



nature of the local current fluctuations around its aver- 



age (Fourier's law) behavior. The functional in (2.13) 
is the same as in the conservative case (that is, with no 
bulk dissipation) , due to the quasi-elasticity of the micro- 
scopic dynamics, which makes the current noise be the 
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FIG. 1. (Color online) Convergence of the space&time- 
integrated dissipation to its ensemble value for many different 
realizations, and sketch of the probability concentration as 
time increases, associated with the large deviation principle, 
eq. |2U5|. 



only relevant one in the hydrodynamic description, see 
discussion in section IV 39 . We focus now on the fluc- 



tuations of the dissipated energy, integrated over space 
and time 



/•1/2 

dt / dx d(x, t) 

J-l/2 
,1/2 

dt / dxR(p(x,t)) > 0, 

(l J-l/2 



(2.14) 



where we have introduced a minus sign for the sake of 
convenience, in order to make d positive. As discussed 
above, this is a fundamental observable to understand 
the statistical physics of driven dissipative media. The 
probability of such a fluctuation P T (d) scales in the long- 
time limit as 



P T (d) ~ cxp [+tL G(d)], G(d) = - max I T [p, j] . 

T p,j 

(2.15) 

This defines a new large deviation principle for d, see Fig. 
[I] such that G(d) is obtained from I T [p, j] via a saddle- 
point calculation for long times (that is, it follows from 
the contraction of the original rate function I T [TO])- The 
optimal fields po(x, t; d), jo(x, t; d) which are the solution 
of the variational problem ( 2.15[ ) must be consistent with 
the prescribed value of the dissipated energy d in (2.14), 
and are also related by the balance equation (2.1), sup- 



plemented with (2.4) and the appropriate boundary con- 
ditions. These optimal fields can be interpreted as the 
ones adopted by the system to sustain a long-time fluc- 
tuation of the space&time-integrated dissipation d. For 
the sake of simplicity, we have not explicitly introduced 
in our notation the parametric dependence of the LDF 
G(d) and the associated optimal profiles on the bound- 
ary temperature T, though this should be borne in mind 
for latter reference. 
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A. The constrained variational problem 



where 



We now assume that these optimal profiles do not de- 
pend on time. In conservative systems, this conjecture 
has been shown [7] to be equivalent to the additivity prin- 
ciple recently introduced to study current fluctuations in 
diffusive media [8 . The validity of this additivity sce- 
nario has been recently confirmed in extensive numerical 
simulations for a broad interval of fluctuations [TTJ Q3] , 
though it may eventually break down for extreme fluctu- 
ations via a dynamic phase transition |30[ 13 Xj . As we will 
see below, the applicability of this generalization of the 
additivity conjecture to dissipative systems is well sup- 
ported by numerical evidence. Under this simplifying hy- 
pothesis, the fluctuating balance equation (2.1 1 reduces 
to 



j'{x) + vy(x) = 0, y(x) = -j'(x)/v, 



(2.16) 



making use of the variable y defined in eq. (2.5a). More- 
over, wc can integrate over time in the definition (2.14) 
of the integrated dissipation d, 



fl/2 

d = v I dx y(x), 

-1/2 



(2.17a) 



or, equivalently, 

,1/2 



d = - 



/ dxj'(x)=j(-l/2)-j(l/2)>0. (2.17b) 

J-X/2 



In this way, by using the additivity hypothesis we can 
eliminate p(x) and write G(d) in terms of only one vari- 
able as 



G(d) = - mmS[j], with S[j] = f dxL{j,j'j") , 
n x > J -i/2 

(2.18a) 



C(j,j'J") = 



2a{-j'/u) 



(2.18b) 



where D is the effective diffusivity d efine d in eq. (2.9|, 
and <t is the mobility, defined in eq. (2.3), both written 
in terms of y = —j'jv. The function C(j,j',j") is a 
generalized Lagrangian with dependence on first and also 
second order derivatives, see the appendix. 

We have to find the optimal current profile jo(x;d), 
that is, the solution of the variational problem (2.18), 



with the constraint that the integrated dissipation d has 
a definite value, as given by (2.17). Therefore we must 



use the Lagrange multiplier procedure [3S[ [33], that is, 
look for an extremum of 



1/2 



(2.19) 



S x [j]=S\j}-\ I dx(j' + d) 

J -1/2 

.1/2 

dxC x (j,j',j") (2.20) 

-1/2 



Cx(j,j',j") = mj'J")-Hf + d), 



(2.21) 



with A being the Lagrange multiplier. The extremum of 
S\ follows from two conditions: (i) SS X — 0, and (ii) 
dS\/d\ = 0. The first condition implies 



dx 2 \ dj" 



± (dLx 
dx \ dj' 



dj 



0, 



(2.22) 



which is the Euler-Lagrange equation for a Lagrangian 
L\ containing second order derivatives (see the ap- 
pendix). Condition (ii) leads to the constraint on the 
integrated dissipation, given by Eq. (2.17). The bound- 



ary conditions for the Euler-Lagrange equation are 

j'{±l/2) = -uR{T), p Ai (±l/2) = 0. (2.23) 

We have introduced the generalized momentum pj con- 
jugate of j, for the new Lagrangian C\, as 



PXj 



dC\ d ( dC), 



dj' dx \ dj" 



(2.24) 



These boundary conditions arise from (1) the values 
of the density at the boundaries, which are prescribed, 
p(±l/2) = T, and (2) the condition SS X = 0, which pro- 
vides the additional needed conditions when there are 
not enough values of the variables fixed at the bound- 
aries (see the appendix, and also [53]). 



B. Mapping the constraint to boundary conditions 

Taking into account the relation between £\ and C, 
Eq. (2.21), the generalized momentum p\j verifies 



p Xj = Pj - A, 



(2.25) 



where pj is the generalized momentum for the Lagrangian 
£, that is, 



Pj 



dC _ d ( dC 
df~ dx \df r 



(2.26) 



Moreover, the Euler-Lagrange equation (2.22) implies 
that 



d 2 ( dC\ d f dC\ dC 



dx 2 \dj" ) dx \dj' ) dj 



= 0, (2.27) 



that is, we also obtain the Euler-Lagrange equation cor- 
responding to the original Lagrangian C. The definition 
of pj in eq. (2.26) guarantees that p'j — dC/dj, as in the 
case of the usual variational problem with a Lagrangian 
including only first-order derivatives. Now, the boundary 
conditions can be written as 



i'(±l/2) = -vR(T), P j(±l/2) = X, 



(2.28) 



() 



which follow from Eqs. (2.231 and (2.25). The above 



result imply that our constrained variational problem 
can be mapped onto a unconstrained variational prob- 
lem with the original Lagrangian C, its associated Euler- 
Lagrange equation (2.27) and the boundary conditions 



(2.28). The unknown value A for the generalized momen- 



tum pj at the boundaries must be determined by impos- 
ing the prescribed value of the integrated dissipation, as 
given by Eq. (2.17), that is, A = X(d). In particular, 
A = is equivalent to imposing no restrictions on the 
integrated dissipation, so that we should recover (as we 
will see later) the average profiles and dissipation in this 
case. In this sense, a non-zero value of A = p., (±1/2) is a 
measure of the departure from the average hydrodynamic 
behaviour. 

On physical grounds, we expect the corresponding op- 
timal density profile to be an even function of x, because 
of the symmetry of our system around the center x = 0. 
In fact, the Euler-Lagrange equation (2.27) admits so- 
lutions with well-defined parity. Since the Lagrangian 
has t he sym metry property C(-jJ', -j") = C(j,j',j"), 
Eq. (2.27) has solutions with j being an odd function 



of x, which implies that y (and therefore p) is an even 
function of x. From now on, we will restrict ourselves 
to these symmetric solutions of the variational problem. 



Thus, Eq. (2.17) reduces to 



2j(-l/2;d) = -2j(l/2;d), 



(2.29) 



so the boundary conditions for the Euler-Lagrange equa- 
tion boil down to 

j'(±l/2; d) = -uR{T), j(-l/2; d) = -j(l/2; d) = d/2 , 

(2.30) 

i.e. much simpler than Eq. (2.28). This follows from pj 
being an even function of x for the solutions with well- 
defined parity we are considering. In this way, symmetry 
considerations lead to the simpler boundary conditions 



(2.30), which in turn allow us to get rid of the Lagrange 
multiplier A. In summary, we have mapped our original 
variational problem with the subsidiary condition that 
the dissipation has a given value to an unconstrained vari- 
ational problem, with the original Lagrangian C(j,j',j") 
and prescribed values of j and j' at the boundaries. 

Once the optimal current profile is obtained, the op- 
timal density profile can be calculated from the balance 
equation (12. 161). Of course, the density profile so ob- 



tained obeys the boundary conditions p(±l/2;d) = T. 
It must be stressed that the Euler-Lagrange equation 
(2.27) with boundary conditions (2.30) gives the correct 



solution to the constrained variational problem when the 
optimal profiles have a well-defined parity. Nevertheless, 
one cannot rule out the existence of symmetry-breaking 
solutions without well-defined parity, since in general a 
variational problem may have multiple solutions |32j . In 
that case, one must solve the more complex variational 
problem comprising the Euler-Lagrange equation (2.27) 



with the boundary conditions (2.28), where the Lagrange 
multiplier A = A(d) is determined by imposing the con- 



strongly supports the validity of symmetric solutions (see 
below) . 

The LDF G(d) depends on d and T through the bound- 
ary conditions (recall that its T-dependence has been 
omitted in notation for simplicity). The derivation of 
the Euler-Lagrange equation ( |2.27 l (see the appendix) 
shows that, for the solutions with well-defined parity, 

5G(d)=p j (l/2)5d + 2p j ,(l/2)v C ^p-5T, (2.31) 

where pj is the generalized momentum conjugate of j, de- 
fined in Eq. (2.26), and pj/ is the generalized momentum 
conjugate of j , 



Pj' 



dC 



(2.32) 



which is an odd function of x for the solutions with well- 
defined parity. In this way Eq. ( |2.31 ) offers a geometric 
interpretation for the values of the generalized momenta 
at the boundaries, as they are directly related to the par- 
tial derivatives of the LDF, 



dG 
~dd 



= Pi(l/2), 



° G = 2,«W). (2.33) 



dT 



dT 



C. A Hamiltonian formulation of the problem 

We will not write the detailed form of the general 
fourth-order differential equation (2.27) for the optimal 



profile j(x;d), since it is not particularly illuminating. 
Instead we now write a set of four equivalent first order 
differential equations arising in the equivalent "Hamilto- 
nian" description. In the following, we sketch the proce- 
dure to introduce the Hamiltonian for a Lagrangian with 
higher-order derivatives [351 [33J , adapted to the present 
case. As the Euler-Lagrange equation is a fourth-order 
differential equation, we should have two canonical coor- 
dinates and their two corresponding canonical momenta. 
The first canonical coordinate is the current j, and we 
choose the second one to be y, which is proportional to 
j' , as given by Eq. (2.16). This choice is suggested by 



the structure of the Lagrangian in eq. (2.18). It is worth 



recalling that the density profile can be directly obtained 
from y by making use of its definition, Eq. (2.5a). Next, 



we introduce the canonical momenta p y and pj conjugate 
to y and j, respectively. The momentum pj has been de- 
fined in Eq. (2.26), and p y is given by 



Py 



dC 
'W' 



(2.34) 



which follows from the definition of py, Eq. (2.32), and 



Eq. (2.16) for y. The Hamiltonian is then introduced in 



the usual way, 



straint (2.17). We note however that numerical evidence 



W = y'py + j'pj -£ = y'py - vypj - C. 



(2.35) 
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After some algebra, we get 

H = \Q{y)pl ~ D-^ViOPv ~ "VPi, (2.36a) 



Q(y) 



= D 2 (y) ' 



(2.36b) 



where we have defined the auxiliary function Q(y), with 
Q(y) > for all y. We have also made use of eq. ( |2 . 16 1, 
of the expression for p y which follows from its definition 
(pl, 



Py = D{v) 
and of the Lagrangian 



j + D(y)y' 
Hv) 



i 



Q(y)p 2 y , 



(2.37) 



(2.38) 



2Di(y) 2 

written in terms of the canonical variables, with the aid 



of eq. (2.37). As usual, H is a function of (y, j,P y ,Pj), 
which satisfy the following set of four "canonical" first- 
order differential equations, 



m 

d P y 



Q(y)p y -D- 1 (y)j, 



cm 

dp, 



= -vy, 



(2.39a) 



(2.39b) 



&H _ dQ(y) P 2 y dD-^y) . 



dy 



dy 2 



dy 



jPy + vpj, (2.39c) 



Pi 



dU 
dj 



D- 1 {y)p y 



(2.39d) 



that are equivalent to the fourth-order Euler-Lagrange 
equation (2.27). Note that, as is usual in physics and in 
order not to clutter our formulae, we have dropped the 
subindex for the optimal profiles, which are now solu- 
tions of the above canonical equations; the same notation 
is used for the canonical variables in the Hamiltonian and 
for the solutions of Hamilton's equations. On the other 
hand, as the Hamiltonian does not depend explicitly on 
x, it is a first integral of the system (2.39): H = const, 
over any of its solutions. This property may be used to 
simplify the integration of the system. 

In general, for a given value of the dissipation d, we 



have to solve the system of equations (2.39) with the 
boundary conditions 

y(±l/2) = R(T), j(-l/2) = -j(l/2) = d/2, (2.40) 



which follow from Eqs. (2.30) and (|2.16|. Again, we 



have to look for solutions of Eq. (2.391) with well-defined 



parity, that is, y and j are even and odd functions of x, 
respectively (and, therefore, p y is odd and pj even). The 
solution of these canonical equations is then inserted into 
the expression of the LDF G(d), which can be written in 
terms of the canonical variables as 



G(d) 



,1/2 1 pl/2 

/ dxC = -- dxQ{y)p 2 y , (2.41) 

J-l/2 1 J-l/2 



«l/2 



by combining eqs. (2.18) and (2.38). In this way, we 



obtain the LDF for an arbitrary value of the integrated 
dissipation d within the Hamiltonian formulation of the 
variational problem. Equation (2.41) shows clearly that 



the most probable (average) profiles correspond to a so- 
lution with p y — for all x, for which G(d) vanishes. 
By substituting p y = in Eqs. ( 2.39c )-( 2.39d ) , we also 
have that pj = for all x. Moreover, Eqs. ( |2.39a ) and 



(2.39b) simplify to eq. (2.8), that is, the average profiles 
are reobtained. Therefore, there is always a solution of 
the canonical equations (2.39) with identically vanishing 



canonical momenta, which corresponds to the average 
solution of the hydrodynamic equation ( |2.8| 4"0] . These 
average hydrodynamic profiles {p a v, jav} lead to the av- 
erage value of the integrated dissipation 



,1/2 

d w = v I dxR(p m ) = v j dxy av (x) 

J-l/2 J-l/2 



,1/2 



1/2 

2 Jav (-l/2). 



(2.42) 



This discussion is consistent with the one below Eq. 



(2.28), which was done within the framework of the 



equivalent Lagrangian description. Fluctuations involve 
non-zero values for the canonical momenta, whose magni- 
tude is then a measure of the departure from the average 
behaviour (d — d av )/d av . 



III. ANALYSIS OF THE LDF IN SOME 
LIMITING CASES 

In the following subsections, we further analyze the 
form of the LDF in certain limits of interest, for which 
some general results can be obtained. First, we focus on 
the behavior of G(d) for small fluctuations around the 
average, where a quadratic shape of the LDF is expected 
(corresponding to gaussian fluctuations). We then ana- 
lyze the limit of weakly-dissipative systems, i/<l, for 
which an adequate perturbative expansion allows us to 
obtain a non-trivial and interesting scaling form for the 
LDF. Finally, we consider the opposite limit of strongly- 
dissipative systems, v 3> 1, for which a different scaling 
for the LDF is found. 



A. Small fluctuations around the average 

As the average behavior corresponds to the particular 
solution of the canonical equations corresponding to van- 
ishing momenta pj = 0, p p = 0, small fluctuations can be 
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thus analyzed by assuming that the canonical momenta 
are small. Let us define the dimensionles parameter 



d — d,. 



(3.1) 



to measure the separation from the average integrated 
dissipation d av . As we have just discussed, the canonical 
momenta vanish for e = 0. We write 

V = y av + eAy, j = j av + eAj, p y = eAp y , pj = eAp }1 

(3.2) 



and linearize Eqs. (2.391 around the average solution, 



that is, we only retain terms linear in e. Then, 
Ay' = Q(y av )A Pv - D-\y av )Aj - j av ,} Um> Ay, 



dy a 



Af = -vAy, 



Ap " =J -^iy^ Ap v 



vAp^ 



Ap\ = D-^y^Apy 



(3.3a) 
(3.3b) 

(3.3c) 

(3.3d) 



The boundary conditions for these equations are 
Ay(±l/2) = 0, Aj(-l/2) = -Aj(l/2) = d w /2. The 
solution of this system of equations must be inserted in 



the large deviation function (2.18). Using the expression 



\2Al\ for the LDF, it follows that 

1/2 
-1/2 



e 2 r 1/2 

G(d) ~ — ;r / dxQ(y av )Ap; 

2 7-1/2 



(3.4) 



for small fluctuations of the dissipation around the av- 



erage. Taking into account (3.1) and the large deviation 



principle (2.15), eq. (3.4) means that the probability of 
such small fluctuations of the integrated dissipation d is 
approximately gaussian, 



Pr(d) 

with A 2 given by 
A 2 = 



cx exp 



-Lt 



(d - rfav) 2 
2d 2 ,A 2 



,1/2 



1-1/2 



dxQ{y m )Ap 2 y 



(3.5) 



(3.6) 



In this way, the gaussian estimation for the standard de- 
viation of the dissipation, by comparing (3.5) to (3.6), is 
given by x = d &v k v /\rL. In order to make a more de- 
tailed study of the LDF, concrete functional dependences 
of the diffusivity D, the mobility a and the dissipation R 
on the density p must be considered. This is done in the 
following sections of the paper, where we will consider 
a broad family of models for which the transport coef- 
ficients can be explicitly obtained. On the other hand, 



it is important to notice that gaussian statistics is only 
expected for small fluctuations around the average dissi- 
pation. In general, the solution of the variational problem 
given by the integration of eq.( 2.39| ), when inserted into 
(2.41), will give rise to non-gaussian statistics (that is, a 
non-quadratic dependence of the LDF) for an arbitrary 
fluctuation of the dissipated energy d. 



B. Weakly-dissipative systems, i/< 1 

We proceed now by analysing the canonical equations 
(2.39) in the limit v <C 1. Unsurprisingly, a regular per- 
turbation expansion in powers of v breaks down, since 



it is not possible to impose the necessary boundary con- 
ditions for the current. This singularity of the elastic 
limit was to be expected on a physical basis, as it is not 
possible to obtain the behavior of weakly dissipative sys- 
tems {y < 1) as a correction around the conservative 
case v = 0, for which p(x) = T and j(x) = 0. Therefore, 
a singular perturbation analysis should be done, looking 
for a suitable rescaling of the variables for v <C 1. Equa- 
tion ( |2.8[ ) for the averages implies that 

y av = R(T) + 0(u), ]av = -uR{T)x + 0{u 2 ), v < 1. 

(3.7) 

The average current vanishes linearly in v in the limit 
v — > + , as expected. Moreo ver, t he ave rage dissipation, 
obtained by combining eqs. (2.42) and (3.7), is given by 



d &v = vR{T) + 0{v 2 ). 



(3.8) 



Therefore, it is sensible to propose the following rescaling 
of variables 



j(x) = vip(x), pj(x) 



U^,(x) 



v 



(3.9) 



which is consistent with the canonical equations (2.39), 
since 



, _ i ., _ i_ an _ du 

v v dpj dTL^ 



m _&h 

dj dip ' 



(3.10a) 
(3.10b) 



with the same Hamiltonian H. In other words, eq. (3.9) 



defines a "canonical transformation" from the pair of 
canonical conjugate variables {j,Pj} to {^,11^}, a trans- 
formation that heals the singular behavior in the v —> + 
limit. The Hamiltonian can be now written as 



n = ^Q(y) P 2 v 



yll^j - vD(y) tpp y 



(3.11) 



in the rescaled variables. Notice that the transformation 
introduced is essential to obtain the correct "dominant 
balance" [42] to the lowest order. In particular, before 
the rescaling, the term proportional to ypj was of the 
order of v and the term proportional to jp y was of the 
order of unity; after the rescaling the orders of magnitude 
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are interchanged, the term proportional to yll^ is of the 
order of unity while the term proportional to ipp y is of the 
order of v. We now start from the zero-th order rescaled 



Hamiltonian by putting v = in eq. (3.11 1, 



1 



= 7 ) Q{y)p y - 



from which we we arrive at 
, dHo 



y 



Q(y)P: 



(3.12) 



(3.13a) 



ijj ■ 



i>> = 



n' = 



= -y, 



= o. 



(3.13b) 
(3.13c) 



(3.13d) 



dp y 

_&Ho 
dy 
dHp 

dH 

' dip 

In order not to clutter our formulas, we do not introduce 
a different notation for the canonical variables, although 
the approximate canonical equations (with Ho) are dif- 
ferent from the exact ones (with H). We have only to 
remember that our results are valid to the lowest order 
in v. The canonical equations ( |3. 13 ) have to be solved 
with the boundary conditions 



y(±l/2) = R(T), 
where 



lK-1/2) = -V>(l/2) = A/2, 

(3.14) 



A= d =R(T)± 

V 0t a . v 



(3.15) 



is assumed to be of the order of unity, that is, d = 0(y) 
or d/d av = 0(1). Thus, our rescaling allows us to obtain 
a solution for the optimal profiles for the density p, by 
inverting the relation y — R(p), and the current j = 
vip, for integrated dissipations d very different from its 
average value d av . 

It is worth noticing that ip is a cyclic variable and its 
conjugate momentum is thus constant, LT^, = IL^o = 
const., see Eq. (3.14); this fact allows us to obtain a 



closed first order differential equation for y(x) in the 
v <C 1 limit. Moreover, by recalling Eq. (2.33), we have 
that 



n 



dG 
9A' 



(3.16) 



which gives the physical interpretation of this first in- 
tegral of the approximate canonical equations: it is the 
partial derivative of the LDF with respect to the rescaled 
dissipation. The Hamiltonian Ho is also constant, since 



it does not depend explicitly on x, and combining (3.12) 
and ( 3.13a| ), 



y'=2Q(y)(H + yIl 4 ,o), y(±l/2) = R(T). (3.17) 

Once this is solved , the rescaled current ip can be ob- 
tained from ( 3.13c[ ) 



V/ = -y, iK-1/2) - -VKV2) = R(T) 



2cL 



(3.18) 



so that the two constants Ho and II^,o will be given in 
terms of the temperature T and d/d av . There are no more 



constants to be adjusted in the solution of eqs. ( 3.17 ) and 



(3.18 ) due to the parity properties of (y, ip): y is and even 
function of x and ip is an odd function of x in the interval 
[—1/2,1/2]. Of course, the average profiles y av (x) and 
j av (x) are reobtained from the canonical equations by 
putting n^, = and p y = therein. Equation (3.12) 
implies that Ho = over the average profiles. 



The simple form of the differential equation (3.17) al- 



lows us to infer some of the properties of the optimal 
profile y(x) associated to a given dissipation fluctua- 
tion in the limit of weakly-dissipative systems. First, 
notice that in general the solution of eq. (3.17) will 



be non-monotonic, exhibiting extrema in the interval 
x E [—5, 5]. Moreover, taking into account that the func- 
tion Q(y) is positive defined, it follows that the profile at 
the extrema will take an unique value 



Vo 



Ho 

n^o 



(3.19) 



Note that II^,o 7^ for d 7^ d av and, moreover, it must 
have a different sign that Ho, i.e. sgn(IL^o) 7^ s gn(iJ ), 
since y{x) > Vx. Therefore, the optimal profile y{x) 
can only have a single extremum (minimum or maxi- 
mum) [41] , that is located at x = because of symmetry 
reasons. By rewriting eq. (3.17) as 



y u = 2Q(y)Ho 1 - 



y_ 

1/0 



(3.20) 



we conclude that the constant Ho and yo~y(x) must have 
the same sign Vx € [—5, 5]- Thus, for Ho > the profile 
y(x) has a single maximum, y{x) > y(±l/2) = R(T) \/x, 
and thus d > c? av . On the other hand, Ho < implies a 
single minimum, y(x) < R(T) Vx and d < d av . All these 
properties are confirmed below for particular examples, 
both analytically and numerically. 

Interestingly, the leading behavior for the LDF can be 
also easily obtained in terms of the first integrals Ho and 
n^ . In fact 



G{d) ~ / dxQ(y)p 2 y = — dx J^- 

2 7-1/2 zJ-i/2 Q(y) 

(3.21) 



and making use of eq. (3.17) 



G(d) ~-(h + n^ J ^ , 



dx y(x) 



OWo+IWM = -Ho(l 



A 

yo 



(3.22) 



Thus, the remaining task consists in writing the con- 
stants Hq and yo (or equivalently Ho and II^o) in terms 
of the integrated dissipation d and the temperature at 
the boundaries T. Once this is done, the LDF follows 



from the simple expression given by eq. (3.22). Further 



more, we may obtain bounds for the profile extremum yo 



10 



by taking into account that G(d) < for d ^ d av : (i) 
for "Ho > we already know that yo is a maximum and 
this implies that y > A, while (ii) for W. < we already 
know that yo corresponds to a minimum and thus yo < A. 
Interestingly, we can also use Eq. (3.22) together with 
Eq. ( 3.16 ) to obtain a simple relation between Jio, G and 
dGj dd, namely 



Ho — G 



dG 



G(d) 



dG{d) 



It must be stressed that this relation only holds in the 
weakly dissipative system limit, in the sense that the 
macroscopic dissipation coefficient is small, v -C 1, and 
the considered integrated dissipation verifies that d = 
0(1/). 

For many systems of interest, the function Q(y) is typ- 
ically a homogeneous function of y of degree 7, 



Q(cy) = c J Q{y), 



(3.24) 



where c is an arbitrary real number, that is, Q(y) oc y 1 . 
This type of dependence is common to many driven dis- 
sipative media, as for instance the general family of mod- 
els that we will study in section IV [UJ ES] or different 
reaction-diffusion systems [43]. However, it should be 
noted that not all systems obey this homogeneity condi- 
tion, e.g. symmetric simple exclusion processes with dis- 
sipative dynamics have a non-homogeneous Q(y) [9l 137j. 
By introducing the scaling 

y(x) = R(T)Y(x), y = R(T)Y , H = R(T) 2 ->H, 

(3.25) 



Eq. (3.17) is transformed into 



Y\xf = 2HQ(Y) (l - ^M) , Y(±l/2) = 1. 

(3.26) 

This is quite a natural transformation: we scale the vari- 
able y with its value R(T) at the boundaries; besides, we 
will be able to find a physically-relevant scaling variable 
for the dissipation LDF. Accordingly, we also introduce 

ijj(x) = R(T)^{x), &'(x) = -Y(x), (3.27a) 

*(-l/2) = ^l/2) = -^- = -^-. (3.27b) 



2R(T) 2d a 



The (even) solution of Eq. (3.26) has the form 



Y = Y(x,n,Y ). 



(3.28) 



As said before, there are no more integration constants 
when solving Eqs. ( 3.26 )-( 3.27), since Y (resp. ^f) is an 
even (resp. odd) function of x. The boundary condition 
is 



Y(x = l/2,M,Y ) = 1, 



(3.29) 



which implies that Yo — Yo(H), the scaled height Yq is 
only a function of the scaled Hamiltonian T-L. Now, taking 



into account that the (odd) solution of Eq. (3.27) has the 
form 



* = y ), 



we have that 



const. (3.23) Therefore, 



V(x = -1/2,H,Y (H)) 



H = H 



2d,, 



(3.30) 



(3.31) 



(3.32) 



that is, % is a function only of the integrated dissipation 
d relative to its average value d av . 

This observation will be used in what follows to find 
a simple scaling form for the dissipation LDF. In fact, 
equation (3.22) can be readily rewritten as 



R(T)^ 2 G(d) 



d/d av 
Y (H) 



H 



(3.33) 



The equation above gives the general scaling of the LDF 
in the limit of weakly-dissipative systems: since both H, 
and Yo are only functions of d/d av , 



V 



7-2 



G(d) = - 



1 



d/d av 
Yo{%) 



(3.34) 



is only a function of d/d av [we have made use of Eq. ( 3.8 ) 
for R(T)]. This is quite a strong result: it means that, for 
each value of 7, all the curves of (d/d av ) J ~ 2 G(d) plotted 
as a function of d/d av fall on a certain "master" curve 
for all values of the dissipation coefficient v, provided 
that v <C 1 so we are dealing with a weakly dissipative 
system. The only hypothesis is that Q(y) must be a ho- 
mogeneous function of y, with an arbitrary degree 7, a 
rather general assumption satisfied in many cases of in- 
terest (see below). That being said, it is important to 
stress that the differential equation (3.26) for Y(x) con- 
tains the function Q{y), that has 7 as a parameter. Thus, 
both (Y,^) and the rhs of eq. (3.34) also contain 7 as 



a parameter; in principle, different physical models with 
different functions Q(y) have different scaling functions. 
The simplest situation appears for 7 = 2, in that case 
Eq. (3.34) predicts that G(d) is only a function of the 
relative dissipation d/d av , with no additional dependence 
on v. 

Finally, it is also interesting to note that the optimal 
profiles y and ip also have simple scaling forms. In fact, 
eq. (3.28), together with Eq. (3.25), implies that 



-Y(x,H,Yo(H)). 



(3.35) 



On the other hand, eq. (3.30), together with Eqs. ( |3.26| ), 
yields that 



#r) = -^ *(x,H,Y (H)). 
v 



(3.36) 
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Therefore, both y(x) and ip(x) multiplied by v/d^ (that 
is, divided by R(T)) plotted as a function of x collapse 
onto a single curve for constant d/d av for each value of 7 
and all possible v -C 1 in the weakly dissipative regime. 
From eqs. (3.35) and (3.36), the optimal profiles for the 



density and the current are readily obtained, since y = 
R(p) and j = ftp. On the other hand, the first integrals 
Wo and n^ follow from Eqs. (3.19) and (3.25), 



Hq — 



' av 
V 



2-7 



1-7 



H 

Y (U)' 
(3.37) 



C. Strongly-dissipative systems, v 3> 1 

We now proceed to analyze the limit v ^> 1, i.e. the 
limit of strongly-dissipative dynamics. Equation (2.11) 
for the average profile implies that y av (x) develops two 
boundary layers of width l v ~ v~ x l 2 close to x = ±1/2, 
where most of the system energy is localized. This is 
reasonable on physical grounds: for v 3> 1 one expects 
that the injection of energy through the boundaries would 
be limited to a small region near them: most of the 
energy has been dissipated before reaching the bulk of 
the system, effectively decoupling the system into two 
almost-independent halves. This picture can be used to 
simplify the integration of the system of canonical equa- 
tions (2.39): we can restrict ourselves to the half interval 
x G [—1/2, 0] and use the boundary conditions 

y(-l/2) = R(T), i(-l/2) = d/2, y'(0) = 0, j(0) = 0, 

(3.38) 

because of the symmetry of the solutions (y even, j odd). 
Now, by introducing the following rescaling (suggested by 
the typical lengthscale l v ~ i/ -1 / 2 ), 



j = y/vijj, X = -Jv [ x + ~ 



P y = VvTly, (3.39) 



we arrive at the equivalent system of equations 
dy 



dX 



= Q(y)n y -D- 1 (y)^, 
dip 

dx = - y > 



(My 

~d~X 



dQ(y)Ul db-\y) 



dy 



dy 



with the boundary conditions 



y(X = 0) = R(T), 



dy 
dX 



(3.40a) 
(3.40b) 

(3.40c) 

(3.40d) 

0, (3.41a) 



4>{x = 0) 



2^ 



, if>(X = y/v/2) = 0. (3.41b) 



Interestingly, v does not appear explicitly in the rescaled 
canonical equations (3.40 1, but only in the boundary con- 
ditions. Therefore, in the limit v — > 00 we have to solve 



(3.40) with the boundary conditions 



y(X = 0) = R(T), ip{X = 0)=d, (3.42a) 



dy 



lim = 0, lim '0 = 0. 

X^oo dX X-yoo 



where we have defined 



2^' 



(3.42b) 



(3.43) 



which is assumed to be of the order of unity. In fact, 
from Eq. (2.17) one gets 



2v 



1/2 



dxy(x) = 



2v^ / 
Jo 



dXy(X), (3.44) 



that is, 



dXy(X), 



(3.45) 



In this strongly-dissipative regime, the canonical equa- 
tions themselves are not simplified, but a physically ap- 
pealing picture emerges: the system decouples in two in- 
dependent boundary shells of width 0(y/u) close to the 
boundaries, where the rescaled variable X = 0(1) (we 
have restricted ourselves to the semi-interval [-1/2,0], the 
solution in [0,1/2] is found by the symmetry arguments 
already used). Moreover, a simple scaling can be derived 



for the LDF G(d), Eq. (2.411, 



G(d) 



1/2 



dxQ(y)p. 



dXQ(y)U 



(3.46) 

where y and II y are the solutions of (3.40) with the 
boundary conditions (3.42). Therefore, both y and H y 



depend on R(T) and d through the boundary conditions 
and G(d) = ^/vF(R(T),d) where F is a certain function. 

A particularly simple situation appear when both the 
mobility a(y) and the diffusivity D(y) are proportional to 
some power of y, so that (i) the function Q(y) is homoge- 
neous, Q(y) oc y 7 , as in Eq. (3.24), and (ii) the effective 
diffusivity D does not depend on y, D(y) = D = const, 
as discussed in Sec. |n] In fact, this is the case for the 
general class of dissipative models analyzed in Sec. |IV| 
The average dissipation for v 1 is 



(3.47) 



2VvDR{T) 



DR(T) 



The canonical equations (3.40) can be analyzed following 



a line of reasoning similar to the one used in the weakly 
dissipative system limit. Since the details are not neces- 
sary for the work presented here, we only give the final 
result for the LDF, that is 



G(d) 



•^F[[R(T)r 2 



(3.48) 
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where J 7 is a certain scaling function. The scaling in Eq. 
(3.48) is more complex than in the weakly dissipative 
system limit. For instance, in the case 7 = 2 we get 

that G(d) — J v/D d/d av ), so the LDF curves, once 

rescaled by (D/v) 1 ^ 2 , collapse for all v ^> 1 when plotted 
as a function of the relative dissipation d/d av . The factor 
y/v in front of the scaling function accounts for the strong 
supression of the fluctuations of the dissipation that takes 
place in strongly-dissipative systems: for a given value of 
the relative dissipation d/d avi the probability of such a 
fluctuation decreases exponentially with y/v. 



pnLTr 



FIG. 2. The model is defined on lattice sites, each one charac- 
terized by an energy pi . The dynamics is stochastic and pro- 
ceeds via random collisions between nearest neighbors where 
part of the pair energy is dissipated to the environment and 
the rest is randomly redistributed within the pair. Such dy- 
namics mimics at the mesoscopic level the evolution of a wide 
class of systems characterized by a nontrivial competition be- 
tween diffusion and dissipation. 



IV. A GENERAL CLASS OF NONLINEAR 
DRIVEN DISSIPATIVE MODELS 

In order to investigate in detail the validity of the gen- 
eral framework presented in previous sections, we now 
introduce a broad class of dissipative lattice models with 
stochastic microscopic dynamics that contain the essen- 
tial ingredients characterizing many dissipative media, 
namely: (i) nonlinear diffusive dynamics, (ii) bulk dissi- 
pation, and (iii) boundary driving. For the sake of sim- 
plicity, we will present them for the one-dimensional (ID) 
case, but the extension to arbitrary dimension is straight- 
forward. 

We thus consider a system defined on a ID lattice with 
N sites. A configuration at a given time step p is given 
by p = {pi,p}, I = 1,...,N, where pi tP > is the en- 
ergy of the l-th site at time p, so the total energyof the 
system at this time is E p = Xw=i Pi,p> see Fig. Kl The 
dynamics is stochastic and sequential, and proceeds via 
collisions between nearest neighbors. In an elementary 
step, a nearest neighbor pair of sites (1,1 + 1) interacts 
with probability 



Pl,p(P) = 



^l,p = Pl, P + Pl+i, P i (4-1) 



where / is a given function of the pair energy S;. p , and L 
is the number of possible pairs. Clearly L ~ N, but the 
particular relation depends on the boundary conditions 
imposed (e.g., L = N + 1 for open boundaries while L — 
N for the periodic case). Once a pair is chosen, a certain 
fraction of its energy, namely (1 — a)£; iP , is dissipated to 
the environment, mimicking the energy drain observed 
in real dissipative media. The remaining energy aXz iP is 
then randomly redistributed between both sites, 



PLp+i = z p aY, 



l,p ■ 



pi+\,p+\ = (1 - z p )aT,i tP , (4.2) 



with z p an homogeneously distributed random number 
in the interval [0, 1]. This microscopic random exchange 
mechanism yields nonlinear diffusion at the mesoscale, 
being an accurate representation of the coarse-grained 
local energy dynamics in many dissipative systems. The 
above dynamics defines the evolution of all bulk pairs, 
I = 1, .... N — I. In addition, and depending on the 



boundary conditions imposed, boundary sites might in- 
teract with thermal baths at both ends, possibly at dif- 
ferent temperatures T L (left) and T R (right). In this case 
the dynamics is 

Pi.p+i = z p ot{ex, p + e L ), Pn, p +i = z p a(e N _ p + e R ), 

(4.3) 

when the first (last) site interacts with its neighboring 
thermal reservoir. Here e^, k = L,R, is an energy 
randomly drawn at each step from the canonical dis- 
tribution at temperature T^, that is, with probability 
prob(efc) = TjT 1 exp(— efc/Tfe) (our unit of temperature 
is fixed by making ks — 1), see Fig. [I] We may consider 
instead an isolated system with periodic boundary condi- 
tions, such that L = N and Eqs. ( |4.1[ ) and (4.2) remain 
valid for I = (I = N) with the substitution p p = pn. p 
(PN+i,p = pi,p)- 

The simplest dynamics corresponds to f(T,i tP ) = 1 in 



Eq. (4.1 ) . In this case all (nearest neighbor) pairs collide 



with equal probability P; p = L~ , independently of their 
energy. This choice (together with a = 1 above) corre- 
sponds to the Kipms-Marchioro-Presutti (KMP) model 
of energy transport [23], which can be considered as a 
coarse-grained description of the physics of a large class 
of quasi-lD real diffusive systems. For instance, it is 
one of the very few instances where Fourier's law can 
be rigorously proved [53]. In addition, the KMP model 
has been used to investigate the validity of the additivity 
principle for current fluctuations [5] and the Gallavotti- 
Cohen fluctuation theorem [5] and its generalization in 
refs. [ITl HP!]. Another simple, but physically relevant, 
choice is /(Ej iP ) = so that F; iP ~ E; iP /(2£ , p ) for a 
large system. A variant of this model has been recently 
used to study compact wave propagation in microscopic 
nonlinear diffusion [38 . In general, the models here in- 
troduced can be regarded as a toy description of dense 
granular gases: particles cannot freely move but may col- 
lide with their nearest neighbors, losing a fraction of the 
pair energy and exchanging the rest thereof randomly. 
The inelasticity parameter can be thus considered as the 
analogue to the restitution coefficient in granular systems 
[26j : energy is conserved in the dynamics only for a = 1, 
while it is continuously dissipated for any < a < 1. 
Thus, in an isolated system (without boundary driving) 
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the energy would decrease monotonically in time. How- 
ever, if energy is injected, for instance via coupling to 
boundary thermal baths as described above, a steady 
state will be eventually reached where energy injection 
and dissipation balance each other. The class of models 
here presented is an optimal candidate to study dissipa- 
tion statistics because: (a) one can obtain explicit pre- 
dictions for the LDF, and (b) its simple dynamical rules 
allow for a detailed numerical study. The chances are 
that our results remain valid for more complex dissipa- 
tive media with similar macroscopic dynamics. 

The hydrodynamic evolution laws for this family of 
models have been recently studied in detail in ref. [55] . 
In the large system size limit, both continuous space and 
time variables can be introduced, as well as the rele- 
vant hydrodynamic fields: energy density p{x,t), current 
j(x,t) and dissipation d(x,t). This mesoscopic descrip- 
tion is expected to be valid not only for this particular 
class of models but for many driven dissipative media 
over a certain hydrodynamic time scale, much larger than 
the one characteristic of the microscopic dynamics. Over 
the fast (microscopic) time scale, the system forgets the 
initial conditions and relaxes to a local equilibrium state 
in which all the properties of the system become func- 
tionals of the hydrodynamic fields. Afterwards, over the 
much slower macroscopic time scale, the system even- 
tually approaches the steady state following certain hy- 
drodynamic law. For the family of models introduced in 
this paper, the time evolution of the energy density can 
be shown to obey a fluctuating balance equation of the 
general form 

d t p(x,t) = -d x j(x,t) + d(x,t) , 

which is just the starting point for the generalization of 
Macroscopic Fluctuation Theory to dissipative systems 
developed in previous sections. The first term in the 
rhs of this equation accounts for the diffusive spreading 
of the energy, and it is also present in the conservative 
case, while the second one gives the rate of energy dissi- 
pation in the bulk. It is important to note here that the 
models microscopic dyna mics must be quasi-elastic (with 



(1 — a) ~ L , see eq. (4.8 1 below) in order to ensure 



that both diffusion and dissipation take place over the 
same time scale in the continuum limit. Using a local 
equilibrium approximation, the current and dissipation 
fields can be expressed as functions of the local energy 
density |39) . In particular, the fluctuating current can be 
written as 



j(x,t) = -D{p)d x p + i, 



where the first term is nothing but Fourier's (equivalently 
Fick's) law with a diffusivity Dip). This transport coeffi- 
cient can be explicitely calculated for the general family 
of models here presented [33], obtaining 



D{p) 



drr 7 f(pr 2 )e r , 



(4.5) 



where /(£) is the function defining the microscopic col- 
lision rate, see eq. (4.1). On the other hand, the second 
term £ in eq. (4.4) is a noise perturbation, white and 
gaussian. These gaussian fluctuations are expected to 
emerge for most situations in the appropriate mesoscopic 
limit as a result of a central limit theorem. Microscopic 
interactions can be highly complicated, but the ensuing 
fluctuations of the slow hydrodynamic fields result from 
the sum of an enormous amount of random events at the 
microscale which give rise to Gaussian statistics at the 
mesoscale. In the present proof of the gaussian 

character of the noise can be given, due to the simplicity 
of the class of models considered [35] ■ The current noise 
amplitude is a{p)/L (that is, the noise strength scales as 
L -1 / 2 ), where <r(p) is often referred to as the mobility in 
the literature. This coefficient can be again explicitely 
computed within the local equilibrium approximation, 



3 



dr r 7 f(pr 2 )e 



(4.6) 



Remarkably, a direct inspection of eqs. (4.5 



reveals a simple relation between mobility anc 
a(p) = 2p 2 D(p) , 



and (4.6) 



diffusivity 
(4.7) 



which is nothing but a general fluctuation-dissipation re- 
lation for the dissipative case. In fact, it is the same one 
as in the conservative case, because of the quasi-elasticity 
of the underlying microscopic (stochastic) dynamics, see 
eq. (4.8). On the other hand, the dissipation field can 
be written as d(x, t) = —i/R(p), where v is a macroscopic 
dissipation coefficient which can be related to the inelas- 
ticity of the underlying microscopic dynamics 



1 



2L 2 



(4.8) 



and R(p) is a new transport coefficient, absent in conser- 
vative systems 



R{p) = p / drr 5 f(pr 2 )e~ 



(4.9) 



Interestingly, R(p) can be related to the diffusivity. By 
differentiating eq. (4.9) with respect to p after a change 
of variables z = T^fp, it is found that 



1 dRip) Rip) 



(4.10) 



6 dp 3p 

In fact, given the diffusivity, this equation can be consid- 
ered as a first order differential equation for Rip). The 
solution thereof, with an appropriate boundary condition 
(normally i?(p = 0) = 0), is equivalent to calculate the 



(4.4) integral in eq. ( |4.9[ ) 



It should be stressed that the dissipation field t) 
has no intrinsic noise, so its observed local fluctuations 
are enslaved to those of the density. This stems from the 
subdominant role of the noise affecting the dissipative 
term: its strength scales as L~ 3 / 2 as a consequence of 
the quasi-elasticity of the microscopic dynamics, see eq. 
(4.8), so it is negligible against the current noise in the 



mesoscopic limit 
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V. LARGE DEVIATIONS OF THE DISSIPATED 
ENERGY 

In what follows we will apply the theory developed in 
sections [TT||III| to the general class of dissipative models 
introduced in section IV More concretely, we will re- 
strict ourselves to the following choice for the collision 
rate function 



f(p) 



r(/3 + 3)' 



(5.1) 



that is, f(p) oc pP , with /? > — 3 but otherwise arbitrary. 
We have introduced the constant 2/T(f3 + 3) |H] for the 
sake of convenience, as it simplifies the expressions of the 
transport coefficients, see below. For 8 = 0, f(p) = 1 and 
all the pairs collide with equal probability, independently 
of their energy value. Thus, the dissipative generalization 
of the KMP model introduced in ref. [25] is recovered. 
For 8 = 1, f(p) = p/3 and the colliding pairs are cho- 
sen with probability proportional to their energy. The 
conservative case has been recently analyzed in ■ 

The transport coefficients for this family of models are 
easily calculated. The dissipation coefficient is readily 
obtained from Eq. (4.9), 



R(p) 



T(p + 3) 



P" 



P+l / drr 5+2/3 e -r 2 =p p+l t (52) 



that gives the rationale behind the choice of the propor- 
tionality constant in Eq. (5.1 1. The diffusivity is calcu- 



lated by substituting Eq. (5.2) into (4.10) 



D(p) 



(i 



(5.3) 



Of course, the same result is obtained with eq. (4.5). 



Finally, the mobility cr(p) follows from the fluctuation- 
dissipation relation ( |4.7[ ), which gives it in terms of the 
diffusivity, 



a(p) = 2p 2 D(p) = 



(3 + 3 



3+1 



(5.4) 



Of course, for 8 = the values of the transport coeffi- 
cients of the dissipative version of the KMP model are 
recovered, D(p) = 1/2, a(p) = p 2 , and R(p) — p [25] . 
Interestingly, the algebraic dependence of transport co- 
efficients with the energy density p appears ubiquitously 
in real systems. One example is granular materials [26] . 
where the density field p may be assimilated to the lo- 
cal granular temperature. For the hard sphere model, 
the average collision rate is proportional to the square 
root of the granular temperature. Thus, this granular 
gas case should correspond to j3 — 1/2, and in fact it is 
found that D(p) oc p 1 / 2 while the dissipative term goes 
as R(p) oc p 3 / 2 . The latter is responsible for the algebraic 
decay with time of the granular temperature (Haff's law; 
oc t~ 2 for large times) observed in the homogeneous case 
when the system is isolated [26] . 



Before going into the details, it is convenient to write 
the explicit form of the auxiliary variable y, related to the 
density p, defined in eq.(2.5a). For our family of models, 



y = P 1+p , 



p = yi 



(5.5) 



where we have made use of eq. (4.9 ). Following the nota- 
tion introduced in section [TT] we also need the "effective" 
diffusivity D{y) 1 given by eq. (2.9), and the mobility 
a 



<j{y), both written in terms of y, 



D(y) = D(p(y)) 



dp _ 3 + /3 
Ty ~ 6(l + /?) : 



<%) = <r(p(y)) 



3 + 8 2+g 

—5— v 1+f> - 



(5.6a) 



(5.6b) 



Hence, as anticipated in section ITT] while the "true" mo- 



bility D(p) depends on p, see eq. (4.5), the "effective" 
diffusivity is constant, D(y) = D. This allows us to cal- 
culate explicitly the average profiles for the density and 
the current, using the linearity on y of eq. (2.11 ) above. 
The steady average solution is thus 



cosh ( s 



y av (j:) = T 



cosh . 



. sinh 



cosh . 



(5.7) 



.(5.8) 



Moreover, the average density is readily obtained by com- 
bining (5.5) and (5.7), 



p av {x) = T 



cosh 



cosh . 



1+3 



(5.9) 



The validity of these hydrodynamic predictions has been 
tested via extensive numerical experiments in ref. [39] . 
Notice that one can define an unique natural lengthscale 
associated to a given v from the hydrodynamic profiles 



above, namely l v — yD/v. This is the lengthscale 
over which profiles vary appreciably, and it decreases 



like 



-1/2 



as v grows. In fact this observation sug- 



gests that, in the limit of strongly-dissipative systems 
v 3> 1, boundary energy layers develop localized around 
the thermal baths, effectively decoupling the system in 
two almost-independent parts, an observation that we 
have already used in Section |III C| to obtain a simple 
scaling relation for the dissipation LDF for strongly dis- 
sipative systems v 1. 

We are interested in the probability of a given fluctua- 
tion of the dissipated e nergy d, integrated over space and 
time, as defined in eq. (2.14). As described in sectionfTTl 



this probability obeys a large deviation principle 
V T (d)~exp [+tL G(d)] , 



(5.10) 
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where the large deviation function G(d) is obtained by 
solving a variational problem for the so-called "optimal" 
profiles {po(x;d),jo(x;d}} which sustain the considered 
fluctuation, after an (additivity) conjecture on the time- 
independence of optimal paths. In fact, making use of 
eq. EM, 



G(d) 



1/2 



-1/2 



dxQ(y)p 2 y , 



(5.11) 



with {y(x) , p y (x )} being the solutions of the canonical 
equations ( |2.39 1 with the appropriate boundary condi- 
tions. For the family of models introduced in eq. (5.1) 



for the collision probability, the Hamiltonian defined in 
eq. (2.36a) is 



H 



1 

~D" 



Q(y)p 2 y - -^3Py - vypj. 



(5.12) 



The auxiliary function Q(y) in the equations above has 
been defined in eq. (2.36b), for the case we are analyzing 



Q(y) = 



D 2 



12(1 + B) 2 2+3 
3 + /3 y ' 



(5.13) 



that is, Q{y) is basically the mobility, since D is constant, 
independent of y. Thus, Q{y) is an homogeneous function 
of degree 7, 



Q(y)=cy J , 7 



2 + /3 
1 + jS' 



(5.14) 



with c= 12(l + /3) 2 /(3 + /3). The parameter 7 varies from 
the value 7 = 2 for the case (3 = 0, which corresponds 
to the dissipative KMP model, to 7 = 1 for the limit 
/3 — > 00. The canonical equations determine the optimal 
profiles, 



/ 9H j 

v = dv- y =Q[y)Pv ~b- 



(5.15a) 




(d-d av )/% 

FIG. 3. (Color online) Probability distribution for the dissi- 
pated energy, integrated over space and a long time r, plotted 
versus the reduced variable (d — d av )/x (of unit variance) for 
many different values of v £ [1CP 3 , 10 3 ], for the case /3 = 0. 
Inset: Semilog plot of the same data. In both cases the line 
is the normal distribution. Gaussian statistics is observed for 
typical fluctuations, but the tails already show signs of asym- 
metry. 



a function of y. Of course, just as in the general case, 
the canonical equations (5.15) have a particular solution 
with vanishing canonical momenta, p y = 0, Pj 



0. 



y = 



d' 



J 



-vy, 



(5.17) 



that is the particularization of eq. (2.8) for our family 
of models. The solution thereof are the average profiles 



(5.7|-(5. 



A. Typical fluctuations and gaussian behavior 



cm 

dp, 



= -vy, 



&H _ dQ(y) p\ 



dy 



Pi 



dy 

dH. _ Py 



vp h 



(5.15b) 



(5.15c) 



(5.15d) 



The boundary conditions for this system of differential 
equations are 

y(±l/2) = T 1+ ^, i(-l/2) = -i(l/2) = ~. (5.16) 

The main differe nce with the general system of canonical 
equations (2.39) comes from D being constant and not 



As described in subsection |III A[ gaussian fluctuations 
are expected for small deviations from the average dissi- 
pation b ehav ior, P T (d) « exp[— Lr{d — d av ) 2 /(2c? 2 v Aj)], 
see eq. (3.5). Following the general theory, small fluc- 



tuations correspond to small canonical momenta, since 
the average behavior is obtained for p y — 0, pj = 0. We 
linearize the canonical equations by introducing the pa- 
rameter e = {d — d av )/d av <C 1. Recalling eq. (2.42), and 
making use of eq. (5.7), the average dissipation is 



f 1/2 I - 

d av = v \ dxy av (x) = T 1+,3 V4i/Dtanh . 

J-1/2 



v 

AD 

Writing the canonical variables as their averages plus a 
linear correction in e, see eq. (3.2), we find at first order 



in e the following set of equations 

Ay' = Q(y av )A Py - -lAj 



(5.19a) 
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A? = -i/Ay, 



(5.19b) 



Ap' = i/Ap, 



ApJ = ^A P y, 



(5.19c) 



(5.19d) 



which particularizes eq. (3.3 1 for our family of mod- 
els. The boundary conditions are Ay (±1/2) = 0, 
Aj(-l/2) = -Aj(l/2) = d av /2. The solution of this 
system must be inserted into eq. (3.6), which gives the 
variance of the gaussian distribution, \ 2 = d av A^/Lr. 
The canonical momentum Ap y is directly obtained by 



integrating eqs. (5.19c) and (5.19d) 



Ap y — K sinh 




(5.20) 



since Ap y must be an odd function of a: as a consequence 
of y being even. The constant K is to be determined with 
the aid of the boundary conditions, but this can only 



be done after solving eqs. (5.19a) and (5.19b), having 



(5.20) into eq. (3.6) gives 



previously inserted (5.20) into them. Substitution of eq. 



At = K 




dxQ{y &v )s\-ah 2 




(5.21) 



In order to evaluate the integral, eq. (5.14) for Q(y) 



must be used, Q(y) oc y 1 , with the parameter 7 being a 
function of /?, 1 < 7 < 2. We now analyze the simplest 
choice (i = 0, that is, 7 = 2, that corresponds to the 
dissipative version of the KMP model introduced in [25] . 
In this case, the calculation is straightforward and yields 



A = 



sinh(2V2i/) - 2\/2v 
4 v / 2^sinh 2 ( v / 2^) 



(5.22) 



Interestingly, A 2 V ~ 1/3 independent of v in the limit of 
weakly-dissipative systems 1/ « 1. This can be under- 
stood as a reminiscence of the scaling of G(d) derived in 
In fact, eq. (3.34) tells us that, for 7 = 2, 



III 



section 

G(d) is just a function of d/d av . In the gaussian ap- 
proximation, this implies the convergence of A 2 to a con- 
stant value in the quasi-elastic limit as v — > + . On the 
other hand, Af, ~ {2y/2v)~ 1 for v ^> 1, which is con- 
sistent with the suppression of dissipation fluctuations 
previously found in the strongly inelastic regime, as ex- 
pressed by the general scaling of the LDF given by Eq. 
(3.48). The same qualitative observations apply to other 



values of j3, though the calculation is more convoluted. 

We have tested the above predictions in standard 
Monte Carlo simulations of the dissipative KMP model 
described in this section, for the particular case /3 = 0. 
Figure [3] shows the probability density function (pdf) for 
the dissipated energy, integrated over the whole system 
and over a long time r for many different values of the 




FIG. 4. (Color online) Measured average dissipation and its 
variance as a function of v for /3 — 0. The solid line corre- 
sponds to the theoretical prediction for d av , eq. (5.18), while 



the dashed line is the gaussian est imati on of the dissipation 



variance parameter, A„, see eq. (5.22). The agreement is 
excellent in all cases. Notice in particular the scaling with 
v of both observables in the weakly- and strongly-dissipative 
system limits. 



macroscopic dissipation coefficient v g [10~ 3 , 10 3 ]. In or- 
der to minimize finite-size effects in the measurements, 
we performed simulations for systems with increasing size 
as v grows, L oc l„, in such a way that the number of 
lattice sites per unit typical length is constant and large 
enough so we are within the hydrodynamic regime. Fur- 
thermore, the integration time r = 0(1) for the contin- 
uous, diffusive, timescale over which the hydrodynamic 
predictions should hold [IS]. Standard Monte Carlo sim- 
ulations do not allow us to sample the tails of the distri- 
bution, but they are useful to study the typical fluctua- 
tions around the average we are interested in here (e.g., 
a regime of 5 standard deviations around the average). 
Figure [3] shows that, when plotted against the reduced 
variable z = (d — d av )/x, the distribution P T (z) follows 
approximately a normal distribution for typical fluctu- 
ations. Moreover, all curves for different v collapse in 
this regime. However, even at this standard simulation 
level, it becomes apparent that the tails of the distribu- 
tion (corresponding to moderate dissipation fluctuations) 
deviate from gaussian behavior, see inset in Fig. [3j show- 
ing asymmetric tails and breaking the collapse to gaus- 
sian behavior observed for small fluctuations. As we will 
show below, the analysis of the dissipation LDF shows 
that the large fluctuations statistics is far from gaussian. 

In order to further check our theory, we have also com- 
pared the measured average dissipation and its variance 
with the analytical results above, as a function of the 
macroscopic dissipation coefficient v, varying in a range 
which covers 6 orders of magnitude. Again, we see in 
Fig. [4]that the agreement is excellent in all cases. In par- 
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ticular, the average dissipation grows as v (resp. v 1 ^ 2 ) 
in the weakly (resp. strongly) dissipative system limit, 
while the variance remains constant for 1/ < 1 but de- 
cays as v~ x l 2 for v 1. Remarkably, the gaussian ap- 
proximation for the variance turns out to be an excel- 
lent estimator of the empirical dissipation variance. For 
Lt 1, large fluctuations of the dissipation are very rare 
and most of the probability concentrates in a region of 
width proportional to (Lt) -1 / 2 around the average value, 
a regime described by the gaussian approximation. 



1. Weakly- dissipative systems, v< 1 

We now focus our attention on the analysis of LDF of 
the integrated dissipation for weakly dissipative systems, 
in which v <C 1. In the general framework developed in 
section [iTl we found a scaling property for G(d), as given 
by eq. (pIL 



7-2 

' nx 1 G(d) = - 
v ' 



1 - 



d/d av 
Y (U) 



H. 



(5.23) 



B. Complete fluctuation spectrum for the 
integrated dissipation 

We now investigate the whole spectrum of fluctua- 
tions (both typical and rare) of the integrated dissi- 
pation. Thus, we need to evaluate the LDF G(d) for 
arbitrary values of d, in general not close to its aver- 
age value c! av , both analytically and numerically. Ex- 
ploring in standard simulations the tails of the dissipa- 
tion distribution associated to the nontrivial structure of 
G{d) is an daunting task, since LDFs involve by defini- 
tion exponentially-unlikely rare events, see eq. (2.12|. 
This has been corroborated in Fig. [3j where the dissi- 
pation distribution has been measured directly but we 
are unable to gather enough statistics in the tails of 
the pdf to obtain clear-cut results in the non-gaussian 
regime. A recent series of works have addressed this is- 
sue, developing an efficient method to measure directly 
LDFs in many particle systems [27H29"] . The method is 
based on a modification of the dynamics so that the rare 
events responsible of the large deviation are no longer 
rare [27], and it has been developed for discrete- [27J 
and continuous-time Markov dynamics |28| . For a re- 
cent review, which also discusses Hamiltonian systems, 
see ref. [29]. The method yields the Legendre-Fenchel 
transform of the dissipation LDF, which is usually de- 
fined as fi(s) = m&Xd[G(d) + sd] [TU] 0B]. In particular, 
if Ucc is the transition rate from configuration C to 
C of the associated stochastic process, the modified dy- 
namics is defined as Ucc(s) = [/c'c ex PNcc), where 
dec is the energy dissipated in the elementary transition 
C -+C '. It can be then shown pTjEZrEH] that the nat- 
ural logarithm of the largest eigenvalue of matrix U(s) 
gives /Lt(s), which in turn can be Legendre-transformed 
back to obtain a Monte Carlo estimate of G(d). The 
method of refs. [27U29] thus provides a way to measure 
/Lt(s) by evolving a large number M of copies or clones 
of the system using the modified dynamics U(s). This 
method is exact in the limit M — ¥ oo, but in practice 
we are able to simulate a large but finite population of 
clones, typically M £ [10 3 ,10 4 ]. This introduces addi- 
tional finite-size effects related to the population of clones 
which must be considered with care, see |13| for further 
discussion along this line. The numerical results for the 
LDF in the following sections have been obtained using 
these advanced Monte Carlo techniques. 



where 7 = (2 + j3)/(l + P),Y (H) is determined by Eq. 
(3.29), and the constant H depends only on the ratio 



d/d av , as given by eq. (3.31 ). 

For the sake of concreteness, let us consider now the 
simplest case j3 = 0, corresponding to the dissipative 
KMP model introduced in [25]. Equation (3.26) for the 
rescaled density profile now reads 

Y\xf = SHY 2 (i-y)> F ( ±1 /2) = 1, (5.24) 
which can be explicitly integrated, with the solution 



Y{x,H) = Y scch 2 (xV2H), Y = cosh 2 \/ j , (5.25) 

where we have already used that Y(x) must be an even 
function of x. The rescaled current profile ^(x) intro- 
duced in (3.27) is 



cosh 



tanh(xV2rO 



(5.26) 



The optimal profiles for the density and the current can 
be now readily written by combining the previous two 



equations with Eqs. (3.35 )-(3.36 1, yielding 



,>{.<■) = TY(.r) = I .sir \ljsech 2 (xVm), (5.27a) 



j(x) = d av ^(x) = -d a 



cosh 



tanh(xV2ft) 



(5.27b) 

in terms of H = H(d). We have taken into account that 
y = p for (3 = 0. Note that the curves p(x)/T — Y{x,'H) 
for different values of v plotted as a function of x only 
depend on the relative dissipation d/d av . Now, eq. (3.31 ) 
implies that 



2*(-l/2) 



(d av - vT), (5.28) 



which gives the constant H implicitly in terms of d/d av . 
Finally, particularizing eq. (5.23) for the case 7 = 2 we 
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FIG. 5. (Color online) Scaling of the dissipation LDF in the 
quasi-elastic limit (y <C 1) for N = 50, T = 1 and varying ft 
for two different values of v, namely v = 0.01 (filled symbols) 
and v = 0.1 (open symbols). The solid lines are the MFT 
predictions in each case. Curves have been shifted vertically 
for convenience, G(d av ) = 0, W, p. For the case /3 = 0, the 
simulation curves are plotted for d < di, with di being the 
inflection point at which G(d) changes convexity in the limit 
1/ < 1 (see the text and also Fig. Rjl. Inset: Comparison 
of the theoretical G(d) for different j3, where it is clear that 
increasing ft favors larger dissipation fluctuations. 



are analyzing (that is, /3 = 0), we obtain 



G(d) = V2-HtanM/— - H. 



(5.29) 



Note that eq. (5.281 for H(d) requires some careful anal- 
ysis. From the general discussion in Secti on |III B| we 

have that % > for d > d av , and thus is a real 

number, while W < for d < d av and y/2H is imaginary. 
The latter case poses no problem for G(d), which is al- 
ways real- valued. In fact, if we write \f% = i\ \H\ we 



arrive at G(d) = -\J2\U\ tan \J ^ + \U\ for H < 0. In 

the limit as d — > 0, we have that T~L —> — 7r 2 /2, and thus 
G(d) — > — oo as expected on a physical basis. 



Equation (5.29) gives a simple scaling form for G(d), 
independent of v, for the linear (/3 = 0) dissipative KMP 
model in the low-dissipation limit i/< 1 25j. As antici- 
pated by (5.23), the curve of G(d) vs the relative dissipa- 
tion d/d av is independent of v in this quasi-elastic regime. 
This scaling is fully confirmed in Fig. [5j in which we plot 
G(d) for different, small values of v € [10^ 2 ,10 _1 ] mea- 
sured in simulations of the dissipative KMP model using 
the advanced Monte Carlo technique described at the be- 
ginning of this subsection. In particular, the agreement 
between theory and simulations is excellent in the broad 
fluctuation regime that we could measure (see below). 
The dissipation LDF is highly skewed with a fast decrease 
for fluctuations d < c? av and no negative branch, so fluctu- 
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FIG. 6. (Color online) Scaling plot of the Legendre transform 
of the dissipation LDF, /i(s) = maxd[G(d) + sd], in the quasi- 
elastic limit v <C 1 for N = 50, T — 1, /? = and two different 
values of v, namely v — 0.01 (circles) and v — 0.1 (triangles). 
The solid line is the MFT prediction, see eq. (5.34 1. Notice 



that n(s) is defined up to a threshold value sj = 0.878458/d av , 
beyond which the Legendre-Fenchel transform diverges. The 
top-right inset shows a zoom around the threshold sj. This 
is related to the existence of an inflection point in G(d) for 
di = 2.27672d av , i.e. a point at which G"(di) = 0, see middle- 
left inset, beyond which the dissipation LDF is non-convex, 
see discussion in main text. 



ation theorem-type relations linking the probabilities of a 
given integrated dissipation d and the inverse event — d do 
not hold [SI [S] • This was of course expected from the lack 
of microreversibility, a basic tenet for the fluctuation the- 
orem to apply [36] . The limit H 3> 1 corresponds to large 
dissipation fluctuations, where G(d) w — i[ln(d/d av )] 2 , 
that is, a very slow decay which shows that such large 
fluctuations are far more probable than expected within 
gaussian statistics (~ — \{d/d &v ) 2 ). In fact, such slow 
decay implies the presence of an inflection point in G(d): 
there is a value di such that G"(di) = 0. The convex- 
ity of G(d) changes at d = dj, G"(d) < for d < dj 
while G"( d) > for d > di 
and (5.29) imply that d]/d e 



Specifically, Eqs. ( |5.28[ ) 
= 2.27672 (see middle- 



left inset in Fig. |6j). The complete measurement of 
non-convex LDFs in computer simulations is a challenge 
which remains unsolved. The reason is that the advanced 
Monte Carlo method described above to directly mea- 
sure large-deviation functions in simulations is based on 
the Legendre-Fenchel transform for the LDF of interest, 
which is not well-behaved in regimes where the LDF is 
non-convex [10] . 

To better understand this issue, recall that the 
Legendre-Fenchel transform of the dissipation LDF can 
be written as 



fi(s) = max[G(d) + sd] = G[d* (s)] + sd* (s) , (5.30) 

d 
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where d* (s) is solution of the equation 
dG(d) _ 



(5.31) 



Note that, mathematically, fi(s) is the Legendre-Fcnchcl 
transform of —G(d), because the Legendre-Fenchel trans- 
form is defined for convex functions |10j . The partial 
derivative of G with respect to d is related to the first in- 
tegral of Hamilton equations H^q, see Eq. (3.16), which 
in turn can be obtained from eqs. (3.371 and (5.25), 
yielding 



n 



■tpO 




Equivalently, 



dG 



dG H l2 H 

•— = — sech \ — . 
3d vT V 2 



(5.32) 



(5.33) 



In this way, making use of Eqs. (5.28), ( 5.29| ) and (5.33), 
the Legendre transform of the dissipation LDF can be 
written as 



H{s) = 2V2£Ttanh' 



-n, 



(5.34) 



in terms of %, which is obtained implicitly as a function 
of s from Eq. (5.33). Note that the scaling of G(d) 
with d/d^v, see Eq. (5.23), implies a similar collapse for 
/i(s) when plotted as a function of sd av . Eq. (5.31) 
has a single solution d*(s) for s < and hence poses 
no problem. On the other hand, due to the existence 
of an inflection point, G'(d) exhibits a minimum at di, 
increasingly smoothly afterward to reach asymptoticaly 
zero in the limit d — > oo, see middle- left inset in Fig. 
[6j Therefore, for s > there exist two solutions d*(s) < 
di < d^i-s) for eq. (5.31 ), but only the first one maximizes 
eq. (5.30). This means that we cannot obtain G(d) by 
inverse Legendre-transforming /Lt(s) for dissipations above 
the inflection point dj — 2.27672 d av . In fact, /i(s) is 
defined up to a critical si, such that sj = 0.87845/d av 
(the slope of —G(d) at the inflection point), beyond which 
fi(s) diverges. This can be seen by noticing the main 
properties of /tt(s), namely 



Oft 

ds 



ds 2 



d 2 G 
3d 2 



(5.35) 



Therefore, \x has a singularity at the value of the slope si 
corresponding to the inflection point dj, where d 2 [i/ds 2 
diverges. The transition to non-convex behavior thus 
implies that we can only measure the statistics of rare 
dissipation fluctuations up to di using the cloning algo- 
rithm [27H29] . Fig. [6] shows a comparison between the 
measured fj,(s) for two different values of v -C 1 and the 
theoretical expectation, up to the critical sj. The agree- 
ment is excellent in all cases, and the collapse of /i(s) 




Po(x;d) 3 



FIG. 7. (Color online) Top: Optimal energy profiles for vary- 
ing d/d a v and f} = 0, measured for v — 10~ 3 (symbols) and 
v = 1CP 2 (dashed lines), and MFT predictions (solid lines). 
Agreement is very good in all cases. Bottom: MFT prediction 
for the optimal density profiles for varying d/d av . 



when plotted against s d av is confirmed. The challenge 
remains to devise computational techniques capable of 
exploring rare-event statistics even in regimes where the 
associated LDF is non-convex. 

We may solve in a similar way the MFT for the in- 
tegrated dissipation for arbitrary values of the exponent 
ft, though mathematical expressions are far more convo- 
luted that in the illustrative case /3 = described above. 
Fig. [5]also shows the dissipation LDF for other exponents 
(3 > 0, as well as the results of numerical experiments in 
these cases. Qualitatively, the results are equivalent to 
those discussed above, with a z^-independent scaling form 
of the LDF in the v <C 1 limit which goes rapidly to zero 
as d — s- and has a relatively fat tail for d ^ d av . This 
tail changes convexity (based on a numerical analysis) 
at a large dissipation di which increases with (a) v for 
fixed (3 (b) (3 for fixed v. For /3=0, we have di/d av ~ 2.8 
for v = 1, while di/d aw > 6 for v = 10. On the other 
hand, for /3 = 1, G"{d) < in the considered region and 
no inflection point therein. For the intermediate value, 
P = 0.5, the positive values of G"(d) are so small that 
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FIG. 8. (Color online) Collapse of the optimal energy profiles 
measured for N = 50 and T — 1 as a function of the relative 
dissipation d/d av for different values of v <C 1, namely v = 
10~ 2 (dotted lines) and v — 10 _1 (dashed lines), for P — 
(top, green), j3 = 0.5 (middle, red) and P = 1 (bottom, blue). 
The larger j3, the less pronounced the central overshoot is for 
d > d av . Solid lines correspond to MFT predictions. 



we have chosen not to eliminate the points behind the 
numerical inflection point, dj/d av ~ 3.2 for 1 and 
di/dav ~ 5.1 for v = 1, although they roughly coincide 
with the values at which the theoretical and the simula- 
tion curves begin to separate. Furthermore, comparison 
with numerical results is excellent in all cases. Interest- 
ingly, see inset in Fig. [5] increasing f3 results in a broader 
dissipation LDF, meaning that large dissipation fluctua- 
tions are enhanced as /? grows away from the linear case 
/3 = 0. 

We have also measured the typical energy profile as- 
sociated to a given dissipation fluctuation for the case 
P = 0, see top panel in Fig. [7j finding also very good 
agreement with the macroscopic fluctuating theory de- 
veloped in this paper. Remarkably, optimal profiles for 
varying v <C 1 also collapse for constant d/d av (all the 
simulations have been done with the same value of the 
energy density at the boundaries T — 1), as predicted by 
eq. (5.27a |. Furthermore, profiles exhibit the x <H- —x 



symmetry conjectured in section [II B| in all cases, with a 
single extremum which can be minimum or maximum de- 
pending on the value of the relative dissipation d/d av , a 
property which was deduced from the general formalism 
in section [ill B| Interestingly, profiles associated to dissi- 
pation fluctuations above the average exhibit an energy 
overshoot in the bulk. This observation suggests that 
the mechanism responsible for large dissipation fluctua- 
tions consists in a continued over-injection of energy from 
the boundary bath, which is transported to and stored 
in the bulk before being dissipated. The same quali- 
tative observations and good agreement between theory 
and simulations is observed for other values of the ex- 
ponent /3 > 0, see Fig. [8j Notice in particular the nice 
collapse of optimal profiles for different values of v <C 1 
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FIG. 9. (Color online) Dissipation LDF for N = 50, T = 1 
and varying /3 = 0, 0.5, 1.0 and v € [10~ 2 , 10]. Curves for /3 = 
0.5 and have been shifted vertically for convenience (recall 
that G(d av ) = W, P), so that P = 1, 0.5 and correspond 
to top (blue), medium (red) and bottom (green). The MFT 
predictions are plotted with lines: solid for P = 1, dashed for 
P — 0.5, and dotted for P = 0. As in Fig. [5] for a fixed v 
increasing P results in larger dissipation fluctuations. 



but equal relative dissipation. An interesting observation 
is that optimal density profiles are less pronounced the 
larger de nonlinearity exponent j3 is, see Fig. [HJ This 
gives a plausible explanation of the widening of G(d) as 
ft increases: for the same value of d/d aw and increasing 
/?, the associated optimal profile is closer to the hydrody- 
namic solution the larger /3 is, and hence this fluctuation 
cost decreases, having a larger associated probability. 



2. Arbitrary dissipation coefficient v 

For arbitrary values of v > 1 no general scaling func- 
tion can be derived in principle for G(d). For each par 
ticular case, the whole variational problem, eqs 
(5.16), must be solved, which is often analytica 



(5.11) 



ly 



tractable. In order to further advance, we resort now 
to a numerical evaluation of the optimal profiles, which 
are used in turn to compute the dissipation LDF. Fig. [9] 
shows the theoretical predictions for G(d) for increasing, 
non-perturbative values of v, together with numerical re- 
sults from simulations, for different values of f3. As for 
the weakly-dissipative system limit previously discussed, 
the agreement between theory and measurements in Fig. 
[9] is quite good. We attribute the observed differences 
between theory and simulation to finite size effects in the 
latter, which are more apparent for large v as compared 
to the weakly-dissipative system limit ^ « 1, compare 
with Fig. [5J see also [39]. Such strong finite-size effects 
are expected since the natural length scale associate d to 

a given v is t v = \l D jv. As follows from Eq. (5.7) 



and the associated discussion, l v decreases as v grows 
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FIG. 10. (Color online) Top: Optimal energy profiles as a 
function of the relative dissipation measured for v — 10, 
N = 50 and T = 1 for the particular case /3 = 0. Thick (green) 
lines correspond to measurements while thin (pink) lines are 
MFT predictions. Bottom: Measured optimal energy profiles 
for v = 10, N — 50 and T = 1, and varying values of the non- 
linearity exponent j3. For a given relative dissipation, energy 
localization around thermal baths decreases as f3 increases. 



so larger system sizes are needed to observe convergence 
to the macroscopic limit. In addition, finite-size effects 
related to the number of clones M used for the sampling 
become an issue in this limit [T3J [3T] • 

In any case, the sharpening of G(d) as v increases 
for any f3 shows that large dissipation fluctuations are 
strongly suppressed in this regime, as was argued for 
v > 1 on quite general grounds in Sec. Ill C I . In this 
strongly-dissipative system limit v ^> 1 the scale t v — > 0, 



and the system decouples effectively into two indepen- 
dent boundary shells. Thus, the energy is concentrated 
around the boundary baths, a picture which agrees again 
with the analysis of Sec. |III C[ This behaviour is evi- 
denced by the optimal energy profiles for a given d mea- 
sured for v = 10, see Fig. [TUJ in contrast to the behavior- 
observed for v -C 1, see Figs. 7]|8 The agreement of the 
observed profiles with MFT predictions is rather good, 
taking into account the non-negligible finite-size effects 
affecting these measurements. Bottom panel in Fig. [T0| 
shows the measured energy profiles as a function of the 
relative dissipation and for different values of the nonlin- 
earity exponent /3. From this figure, it is clear that for 
a given relative dissipation, energy localization around 
thermal baths decreases as (3 increases. This suggests 
again that, as in the v <C 1 limit, the probability of a 
fixed relative dissipation fluctuation d/d av , increases as 
ft grows, giving rise to a broadening of G(d) with /?. 



VI. SUMMARY AND CONCLUSIONS 

In this paper we have developed a general theoreti- 
cal framework for calculating the probability of large de- 
viations for the dissipated energy in a general class of 
nonlinear driven diffusive systems with dissipation. Our 
starting point is a mesoscopic fluctuating hydrodynamic 
theory for the energy density in terms if a few slow hy- 
drodynamic fields, that is, a fluctuating reaction-diffusion 
equation with a drift term compatible with Fourier's law 
and a sink term which can be written in terms of the 
local energy density. The validity of this hydrodynamic 
description can be demonstrated for a large family of 
stochastic microscopic models [35], but it is expected to 
describe the coarse-grained physics of many real systems 
sharing the same main ingredients, namely: (i) nolinear 
diffusive dynamics, (ii) bulk dissipation, and (iii) bound- 
ary driving. From this fluctuating hydrodynamic descrip- 
tion, and using a standard path integral formulation of 
the problem, we can write the probability of a path in 
mesoscopic phase space, that is, the space spanned by 
the slow hydrodynamic fields. Interestingly, the action 
associated to this path, from which large-deviation func- 
tions for macroscopic observables can be derived, has the 
same simple form as in non-dissipative systems. This is 
a consequence of the quasi-elasticity of microscopic dy- 
namics, required in order to have a nontrivial competition 
between diffusion and dissipation at the mesoscale |39j . 

We use the derived action functional to investigate the 
large deviation function of the dissipated energy. The 
energy dissipated in a non-conserving diffusive system 
is, together with the energy current, the relevant macro- 
scopic observable characterizing nonequilibrium behav- 
ior. A simple and powerful additivity conjecture simpli- 
fies the resulting variational problem for the dissipation 
LDF, from which we arrive at Euler-Lagrange equations 
for the optimal density and current fields that sustain 
an arbitrary dissipation fluctuation. A Hamiltonian re- 
formulation of this variational problem greatly simplifies 
the calculations, allowing us to analyze the general the- 
ory in certain interesting limits. A perturbative solution 
thereof shows that the probability distribution of small 
(that is, typical) fluctuations of the dissipated energy is 
always gaussian, as expected from the central limit the- 
orem. Moreover, a general expression for the variance of 
the distribution in the gaussian approximation has been 
derived which compares nicely with numerical results. 
On the other hand, strong separation from the gaussian 
behavior is expected for large dissipation fluctuations, 
with a distribution which shows no negative branch, thus 
violating the Gallavotti-Cohen fluctuation theorem as ex- 
pected from the irreversibility of the dynamics. Further- 
more, the dissipation LDF exhibits simple and general 
scaling forms in the weakly- and strongly-dissipative sys- 
tem limits, which can be analyzed in general without 
knowing the explicit solution of the canonical equations. 

We apply our results to a general class of diffusive 
lattice models for which dissipation, nonlinear diffusion 
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and driving are the key ingredients. The theoretical pre- 
dictions, which can be explicitely worked out in certain 
cases, are compared to extensive numerical simulations of 
the microscopic models (which cover both typical fluctu- 
ations and rare events) , and excellent agreement is found 
in all cases. In particular, the simple scaling for the dissi- 
pation large-deviation function in the weakly-dissipative 
system limit is fully confirmed for different values of the 
nonlinearity exponent /3, exhibiting non-convex behavior 
for large enough fluctuations. Interestingly, in this limit 
energy profiles associated to large dissipation fluc- 
tuations exhibit an overshoot in the bulk resulting from 
an excess energy injection from boundary baths. On 
the other hand, in the strongly-dissipative system limit 
v ^s> 1 the typical lengthscale goes to zero and the system 
decouples into two almost-independent boundary shells, 
giving rise to a different scaling form for the LDF and a 
strong suppression of the dissipation fluctuations in this 
regime. 

Recently, a similar hydrodynamic theory has been de- 
veloped to study large fluctuations in a particular class 
of driven dissipative media [37], but its predictions do 
not compare well with numerical results of the dissipa- 
tive lattice models here studied. The reason is that Ref. 
|37j studies systems with two competing dynamics, one 
conservative and another nonconservative, thus resulting 
in independent fluctuations for the density and dissipa- 
tion fields. In our theory, as is the case in many driven 
dissipative systems, dissipation is linked to the collision 
process, and hence dissipation fluctuations are enslaved 
to density profile deviations. In fact, both theories coin- 
cide in the limit where the optimal dissipation profile is 
given in terms of the optimal density field. 

In summary, our results show that a suitable general- 
ization of macroscopic fluctuation theory [7 is capable 
of describing in detail the fluctuating behavior of gen- 
eral nonlinear driven dissipative media. In this scheme, 
the dissipation LDF follows from a variational problem 
whose solution also gives the optimal profiles that the 
system has to sustain to achieve the considered fluctua- 
tion. The proposed framework is very general, as MFT 
is based only on (a) the knowledge of the conservation 
laws governing a system, which allow to write down the 
balance equations for the fluctuating fields, and (b) a few 
transport coefficients appearing in these fluctuating bal- 
ance equations. This opens the door to further general 
results in the nonequilibrium statistical physics of dissi- 
pative media. In particular, it would be interesting to ex- 



plore the existence of phase transitions and spontaneous 
symmetry breaking at the fluctuating level therein, in a 
way similar to the phenomenon reported in conservative 
systems EI] • Moreover, as the relevant magnitudes 
characterizing nonequilibrium behavior in dissipative sys- 
tems are both the dissipated energy and the current, it 
would be worth analyzing the joint fluctuations of these 
two observables within the MFT approach. 
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Appendix: Variational problem with a Lagrangian 
including second-order derivatives 

Let us analyze a variational problem in which the "ac- 
tion" is defined as the integral of a "Lagrangian" with 
second order derivatives, that is 

S\j]= r dxC(j,j',j"). (A.l) 

The action S[j] is a functional of the profile j(x) in the 
fixed interval x\ < x < xi. The variational problem 
arises when one looks for the "optimal" profile j(x) for 
which the functional S[j] is a extremum. For the sake of 
concreteness, let us consider a problem similar to the one 
analyzed in this paper: we are interested in calculating 
G defined as 

G = -xmn S\j}. (A.2) 

j 0*0 

Then, we consider the variation 5S of the functional when 
a given profile j{x) is slightly changed to j(x) + Sj(x), 

Now, we take into account that 

5f = ~5j, Sj" = ^5j (A.4) 
ax dx z 

in order to integrate by parts (i) once the term propor- 
tional to Sj' (ii) twice the term proportional to Sj". We 
arrive thus at 
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where [f\%\ — f{x2) — f{%i)- By analogy with the case of the usual Lagrangian with only first-order derivatives, we 
introduce the generalized momenta as 



Pj 



d£ _d_ ( dC 
dj' dx \df 



Pj' = 



dC 

W 



In this way, the boundary term has the usual form and Eq. (A. 5 1 can be rewritten as 
SS = \pj8j +Pj'Sj']l 2 i + I dx 



dC _ d /dC\ d 2 (dZ 
dj dx \dj' J dx 2 \dj" 



(A.6) 



(A.7) 



The extremum condition is SS = 0. If the values of j 
and j' are prescribed at the boundaries, both Sj and Sj' 
vanish at x\ t 2 and the boundary term vanishes. Then, 
as Sj is arbitrary for x\ < x < X2, the "optimal" pro- 
file solution of the variational problem verifies the Euler- 
Lagrange equation 



dx 2 



dC 

w 



d_ fdC_ 
dx \dj' 



dC 

dj 



0, 



which is a fourth-order differential equation 
can be written as dpj/dx 



(A.8) 

Interest- 
ingly, Eq. (IA.8I) can be written as dpj/dx = dC/dj that 
is formally identical to the usual Eulcr- Lagrange equation 
for Lagrangians with only first-order derivatives. The 
boundary conditions for the Euler-Lagrange equation are 
the prescribed values of j and j' at the boundaries; 
four conditions for the fourth-order differential equation. 
However, in physical problems there are sometimes less 
prescribed quantities at the boundaries than necessary. 
In that pointed out by Lanczos [33], the ex- 

tremum condition SS = provides the "missing" bound- 
ary conditions. For instance, if we only have fixed values 
of j' at the boundaries (as in the LDF problem we have 
dealt with in the main text), Sj'(x\) — Sj'[x 2 ) = but 
Sj(xi) and Sj(x 2 ) are free parameters. Equation (A.7) 



still implies the Euler-Lagrange equation but also that 



pj(xi) =Pj(x 2 ) = 0. 



The generalized momentum conjugate of the variable 
that is not fixed at the boundary must vanish: the so- 
lution of the variational problem verifies then the Euler- 
Lagrange equation (A.8) with the prescribed values of j' 



at the boundaries and the "extra" conditions provided 
by Eq. (A.9). In this way, we obtain the four condi- 



tions needed to determine completely the solution of the 
Euler-Lagrange equation. 



The function G defined in Eq. (A. 2) depends on the 



values of j and j' at the boundaries. Making use of Eq. 
(A.7), and taking into account that the optimal profile 
j(x) verifies the Euler-Lagrange equation, we get 



$G = ~Pj, 2 Sj 2 -Pj',2Sj' 2 +Pj,iSji 
We have introduced the notation pj 
Sj(xi), 



■Pj 



1,2 and so on. Equation (A. 10) implies that 



Pj(Xi 



(A.10) 

Sji = 



Pj,2 



8G 



Pj',2 



dG 



Pj,i 



8G 
dji ' 



Pj',1 



dG_ 

~~ dj',' 
(ATI) 

Equation ( 2.33 1 of the main paper is the particularization 
of this result for the case (i) X2 = —x% = 1/2, (ii) solu- 
tions of the Euler-Lagrange equation with well-defined 
parity, in which pj 2 — Pj,i, Pj',2 — —Pj',i> and (iii) the 



(A.9) boundary conditions of Eq. (2.30) 
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